# A Study of Coupling Parameter Estimation Implemented by 4D-Var and EnKF with a Simple Coupled System

Authors: Guijun Han; Xinrong Wu; Shaoqing Zhang; Zhengyu Liu; Ionel Michael Navon; Wei Li
Journal: Advances in Meteorology
 (2015)

---

## Abstract

Coupling parameter estimation (CPE) that uses observations to estimate the parameters in a coupled model through error covariance between variables residing in different media may increase the consistency of estimated parameters in an air-sea coupled system. However, it is very challenging to accurately evaluate the error covariance between such variables due to the different characteristic time scales at which flows vary in different media. With a simple Lorenz-atmosphere and slab ocean coupled system that characterizes the interaction of two-timescale media in a coupled “climate” system, this study explores feasibility of the CPE with four-dimensional variational analysis and ensemble Kalman filter within a perfect observing system simulation experiment framework. It is found that both algorithms can improve the representation of air-sea coupling processes through CPE compared to state estimation only. These simple model studies provide some insights when parameter estimation is implemented with a coupled general circulation model for improving climate estimation and prediction initialization.


---

## Body

## 1. Introduction

Due to its potential to reduce initial shocks between different media in a coupled climate system, coupled data assimilation (CDA) that uses coupled model dynamics to extract observational information in one or more media is emerging as an important topic in the climate community ([1–6] the related discussion at the Sixth World Meteorological Organization Data Assimilation Symposium (http://das6.cscamm.umd.edu/)). Based on Bayes’ rule, two main data assimilation algorithms, four-dimensional variational analysis (4D-Var) [7–9], and ensemble Kalman filter (EnKF) [10] have been used to develop CDA systems [3, 4]. Many efforts have been made to compare the performances of 4D-Var and EnKF either under CDA or not [11–16]. It has been found that 4D-Var and EnKF have a comparable performance when the former uses an appropriate minimization time window (MTW) and the latter adopts a suitable variance inflation scheme [17]. While 4D-Var requires a shorter spin-up time in the weather timescale (like days), EnKF can produce a better forecast skill at a long lead time [18–20] in the season timescale (like months).To reduce model errors and improve coupled model predictability, coupled model parameter estimation (also referred to as parameter optimization in the literature) has been introduced into CDA [21–27]. However, coupling parameter estimation (CPE) that uses observations in one medium to estimate the model parameters in other media has not yet been studied in 4D-Var and EnKF. In a climate system that has multiple characteristic time scales, model parameters have remarkable impacts on all model variables. Thus, it is expected that CPE may further enhance the consistency of estimated parameters and model states in a coupled system. Although data assimilation in a multiple space and time scale system has been explored [2, 28], the CPE in such a system has not been fully investigated. Both 4D-Var and EnKF can be used to implement CPE. It is well known that although they originated from the same information estimation theory (Bayes’ rule), different numerical implementations make different performances of 4D-Var and EnKF data assimilation [17]. On one hand, by minimizing a cost function that measures the distance between observations and model states within a specific MTW, 4D-Var CPE seeks a posterior maximum likelihood solution of model parameters in different media in terms of the best fitting of modeling trajectory to observations. On the other hand, EnKF CPE uses flow-dependent coupling error covariance (i.e., error covariance between a model variable in the observational medium and a parameter in another medium) to project observational information onto the parameter being estimated, thereby implementing CPE in a sequential manner.Since evaluating the error covariance of variables residing in different media with a finite ensemble is difficult [6] and a specific MTW is difficult to capture multiple time scales in minimization, CPE is therefore challenging in both EnKF and 4D-Var. A fundamental issue is can 4D-Var or EnKF CPE improve representation of the model air-sea coupling process? Based on a conceptual coupled model [25] which couples a chaotic atmosphere [29] with a slab ocean, we set up an observing system simulation experiment (OSSE) [30] which takes a Nature Run as the truth of model states to answer this question. Within this framework, in both 4D-Var CPE and EnKF CPE, “observations” drawn from a “truth” model that uses the default parameter values are assimilated into the assimilation model which uses erroneously set parameter values for optimizing the coupling parameters.The paper is organized as follows. Section2 describes the methodology, starting from introducing the simple coupled model, followed by the OSSE setup as well as implementation of the 4D-Var CPE and EnKF CPE. This section also presents the scale analysis of perturbed terms in the model equations and the sensitivity study of model parameters. Sections 3 and 4 examine the results of 4D-Var CPE and EnKF CPE, respectively. Different performances of EnKF CPE and 4D-Var CPE are investigated in Section 5. Summary and discussions are given in Section 6.


## 2. Methodology

### 2.1. The Simple Coupled Model

To study coupling parameter estimation (CPE) implemented by 4D-Var and EnKF so as to answer the question (i.e., can 4D-Var or EnKF CPE improve representation of the model air-sea coupling process?) posed in Section1, we employ a simple “climate” model [25] which consists of a chaotic atmosphere model [29] coupled to a slab ocean model to simulate the interaction of the fast atmosphere and slow ocean:(1)
x
˙
1
=
-
σ
x
1
+
σ
x
2
,
x
˙
2
=
-
x
1
x
3
+
1
+
c
1
w
κ
x
1
-
x
2
,
x
˙
3
=
x
1
x
2
-
b
x
3
,
O
m
w
˙
=
c
2
x
2
-
O
d
w
+
S
m
+
S
s
c
o
s
2
π
t
S
p
d
,where an overdot denotes time tendency; x
1, x
2, and x
3 (hereafter denoted by x
1,2
,
3 if they are presented together) are the high-frequency variables of the atmosphere and w represents the slab ocean. (O
m, O
d) define the oceanic time scale, where O
m is the ratio of heat capacity between the slab ocean and the chaotic atmosphere while O
d is the damping coefficient of the temperature w (°C) of the slab ocean. O
m is a dimensionless parameter. The unit of O
d is 1
/
Δ
t, where Δ
t represents the dimensionless time step. The “atmosphere” and the “ocean” interact with each other through the coupling coefficients c
1 and c
2. Here, the unit of c
1 is (
C
°
)
-
1 while c
2 is a dimensionless parameter. The external solar forcing is represented by S
m
+
S
s
c
o
s
⁡
(
2
π
t
/
S
p
d
), where S
m represents the annual mean temperature (°C) of the slab ocean and S
s (in °C) and S
p
d (in Δ
t) represent the amplitude and period of the model seasonal cycle.While the detailed description of the model construction can be found in [25], here, we only comment on the setting of model parameter values. The “atmospheric” parameters σ, κ, and b take their standard values of 9.95, 29, and 8/3. For maintaining the chaotic nature of the “atmosphere” and the stability of the system, the values of c
1 and c
2 are chosen as 0.1 (
C
°
)
-
1 and 1. O
m and O
d are set to 10 and 1, which defines the oceanic time scale as 10 times of the atmospheric time scale. The parameters S
m and S
s are set as 10°C and 1°C, respectively. The S
p
d is chosen as 10 Time Units (TU, 1 TU = 100 time steps = 100Δ
t) so that the period of the forcing is comparable with the oceanic time scale, defining the time scale of the model seasonal cycle. Given the value of S
p
d, the model calendar year is defined as 10 TUs. If we assume that one year has 360 model days, Δ
t is equivalent to 0.36 model days.Using a leap-frog time stepping scheme with a Robert-Asselin time filter [31, 32] and starting from the initial conditions (
x
1
,
x
2
,
x
3
,
w
)
=
(
0
,
1
,
0
,
0
) (note that since the leap-frog time stepping scheme is applied, initial conditions at both time 0 and time 1 are set to (
0,1
,
0,0
).), the model is first freely run for 104 TUs by setting the parameters (σ, κ, b, c
1, c
2, O
m, O
d, S
m, S
s, S
p
d, and γ), where γ is the time filtering coefficient, to be the default values as prescribed in Table 1.Table 1
Default values and values used in OSSE of the model parameters.
Parameter
Default value
Value used in OSSE
Unit
σ
9.95
9.95
—
κ
29
29
—
b
8/3
8/3
—
c
1
0.1
0.11
(°C)−1
c
2
1
1.1
—
O
m
10
10
—
O
d
1
1
1/Δ
t
S
m
10
10
°C
S
s
1
1
°C
S
p
d
10
10
Δ
t
γ
0.25
0.25
—


### 2.2. OSSE Setup

OSSE that takes a Nature Run as the “truth” is an effective way to preliminarily study the validity of data assimilation algorithm and/or assess the observing system impact [30]. In this study, we setup an OSSE to study the EnKF CPE and the 4D-Var CPE. For simplicity, we assume that the model error only comes from the errors of the parameters to be estimated (see Table 1). The coupled model with the default parameter values produces the true solution for the parameter estimation problem. The “observations” are the samples of the “truth” model states after the spin-up of 104 TUs. A random Gaussian noise is superimposed on the “truth” values every 5 time steps for x
1,2
,
3 and every 20 time steps for w, respectively, to produce the corresponding observations. Note that, to simulate the real world sampling frequency, the atmospheric observations are taken more frequently (4 times) than the oceanic observations. The standard deviations of observational errors are 2 for x
1,2
,
3 and 0.2 for w, respectively.The parameters to be estimated arec
1 and c
2 that are the coupling parameters between the ocean and the atmosphere. Through examining the estimation of c
1 and c
2, we can directly study the impact of coupling parameter estimation on the representation of coupling processes. First guess of the parameters c
1 and c
2 is 0.11 (the “truth” is 0.1) and 1.1 (the “truth” is 1.0), respectively. From the initial condition (
x
1
,
x
2
,
x
3
,
w
)
=
(
0
,
1
,
0
,
0
), the assimilation model with the first guess of c
1 and c
2 is spun up for 104 TUs to produce the assimilation initial conditions of x
1,2
,
3 and w (denote as Π). Continuous integrations from these initial values of the assimilation model serve as a control run without any observational constraint. Therefore, this is a free model run, that is, without data assimilation.To study the 4D-Var CPE and EnKF CPE, we design four experiments for each algorithm:(
1
) EXP-1: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean while keeping c
1 and c
2 constant; (
2
) EXP-2: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean as well as c
2 while keeping c
1 constant; (
3
) EXP-3: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean as well as c
1 while keeping c
2 constant; (
4
) EXP-4: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean, c
1 and c
2. Table 2 gives the details of assimilation schemes of these four experiments. The data assimilation period is set to be 100 TUs. Parameter estimation is activated after 10 TUs of state estimation which reaches quasi-equilibrium so that the state-parameter covariance used for parameter estimation is signal-dominated [25]. In addition, the leap-frog time stepping requires two time-level adjustments [33] (i.e., observations at time t are used to adjust the model states at time t and time t
-
1). And also, to investigate the sensitivity of the 4D-Var CPE and the EnKF CPE with respect to the observational interval, four pairs of (atmosphere and ocean) observational intervals: (0.05, 0.2), (0.1, 0.4), (0.2, 0.8), and (0.5, 2.0) TUs are examined. Finally, the results of the last 50 TUs are used for CPE evaluation. Note that, given the low-order model and the scalar nature of parameters, the absolute error of the estimated parameter or the air-sea coupling processes which refer to the truth is the main criterion of the CPE evaluation.Table 2
The instantaneous state and parameter estimations in four data assimilation experiments.
Experiment
Observations
Adjusted variables and parameters
EXP-1
x
1,2
,
3
o, w
o
x
1,2,3, w
EXP-2
x
1,2
,
3
o, w
o
x
1,2,3, w, c
2
EXP-3
x
1,2
,
3
o, w
o
x
1,2,3, w, c
1
EXP-4
x
1,2
,
3
o, w
o
x
1,2,3, w, c
1, c
2We see from time series ofx
2 (Figure 1(a)) and w (Figure 1(b)) in the truth, model control, and “observations” that the erroneous parameters can lead the model states of the free run to depart far away from the truth and observations. According to the governing equation (Equation (1)) of the simple coupled model, the Lorenz-atmosphere is forced by the term c
1
w from the slab ocean while the slab ocean couples with the atmosphere via the term c
2
x
2. Thus, we examine the c
1
w and c
2
x
2 errors to study the representation of the air-sea coupling processes in this study.Figure 1
Time series of (a)x
2 for the first 10 TUs and (b) w for 100 TUs of the assimilation period derived from the truth (red), model control run (blue), and observations (black solid circle), respectively.
(a)
(b)


### 2.3. Scale Analysis of the Perturbed Terms in the Model Equations

Within the OSSE framework, here, we make a scale analysis of the size of the perturbed terms in the model equations. Since the model error is assumed to arise from the uncertainties ofc
2 and c
1, in the ocean, the impact of an erroneous parameter c
2 and an erroneous atmosphere x
2 is contained in the term c
2
x
2. The error is δ
c
2
x
2
=
x
2
δ
c
2
+
c
2
δ
x
2. The characteristic scales of δ
c
2 and δ
x
2 are 0.1 and 10, respectively. Thus, given the default value (i.e., 1.0) of c
2 and the characteristic scale (i.e., 10) of x
2, the scale [i.e., the term x
2
δ
c
2 which is about 1.0] of correction of the oceanic parameter c
2 would be one order smaller than the scale [i.e., the term c
2
δ
x
2 which is about 10] of correction of the atmospheric state x
2.Similarly, in the atmosphere, the impact of an erroneous parameterc
1 and an erroneous ocean w is contained in the term c
1
w. The error is δ
c
1
w
=
w
δ
c
1
+
c
1
δ
w. The characteristic scales of δ
c
1 and δ
w are 0.01 and 1.0, respectively. Thus, given the default value (i.e., 0.1) of c
1 and the characteristic scale (i.e., 10) of w, the scale [i.e., the term w
δ
c
1 which is about 0.1] of correction of the atmospheric parameter c
1 would be the same as the scale [i.e., the term c
1
δ
w which is about 0.1] of correction of the oceanic state w.To summarize, the error of the oceanic state is one order smaller than the error of the atmospheric state. Thus, estimating the oceanic parameterc
2 will be more difficult than estimating the atmospheric parameter c
1.


### 2.4. Model Sensitivities with respect to Parameters

We perform a sensitivity study for seven empirical parameters (σ, κ, b, O
m, O
d, c
1, and c
2) using the same strategy as [26]. Each examined parameter is perturbed through adding a Gaussian noise to its default value (a percentage of the default value as the standard deviation of the noise) while other parameters remain unperturbed. With an ensemble size of 20, the model is freely run for 20 TUs. It is found that, as long as the perturbed parameter values can maintain the stochastic nature of the “atmosphere,” the examined ensemble spread does not have much dependence on the perturbation amplitude (i.e., the standard deviation of the above Gaussian noise) but the spin-up time is a little longer for a small perturbation. For a 5% perturbation (i.e., the standard deviation of the Gaussian noise is 5% of the default value), Figure 2 shows the time series of the ensemble spreads of x
2 (Figure 2(a)) and w (Figure 2(b)), where the black, blue, red, green, yellow, dashed-black, and pink curves, respectively, represent the results of b, c
1, c
2, κ, O
d, O
m, and σ. Here, the ensemble spreads of x
2 and w have been normalized by their own climatological standard deviations (14.5 for x
2 and 1.5 for w). We can see that different variables have different response times (defined as the spin-up period of ensemble spread) with regard to various parameters. For the four atmospheric parameters, the response times of x
2 and w are about 3 TUs and 6 TUs, respectively. For the three oceanic parameters, the response times of x
2 and w are about 5 TUs and 10 TUs, respectively. Therefore, if 4D-Var is used to perform CPE of the oceanic parameters, the typical lengths of MTWs for x
2 and w could, respectively, be 5 TUs and 10 TUs. For example, for x
2, an MTW less than 5 TUs is regarded as a short window while an MTW greater than 5 TUs is regarded as a long window.Figure 2
Time series of sensitivities of (a)x
2 and (b) w with respect to σ (pink), κ (green), b (black), O
m (dashed-black), O
d (yellow), c
1 (blue), and c
2 (red). Note that here the sensitivities of x
2 and w are computed as the ensemble spreads normalized by their climatological standard deviations (i.e., 14.5 and 1.5). Parameters σ, κ, b, O
m, and c
2 are dimensionless while the units of O
d and c
1 are 1
/
Δ
t and (°C)−1, where Δ
t is the dimensionless time step.
(a)
(b)From Figure2, we can see that, among four atmospheric parameters (σ, κ, b and c
1), both x
2 and w have the largest sensitivities with respect to κ (green) and c
1 (blue) during the spin-up period. Similarly, for the three oceanic parameters (O
d, O
m, and c
2), x
2 and w are most sensitive to the coupling parameter c
2 (red). Combining the question presented in Section 1, two coupling parameters c
1 and c
2 are chosen as the parameters to be estimated in the following CPE experiments.


### 2.5. Implementation of 4D-Var CPE

Consistent with Table2, the formulas of the cost functions for four experiments are(2)
J
E
X
P
-
1
X
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
R
-
1
H
X
t
i
-
Y
o
t
i
,
(3)
J
E
X
P
-
2
X
,
c
2
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
V
c
2
c
2
-
c
2
b
2
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,
(4)
J
E
X
P
-
3
X
,
c
1
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
V
c
1
c
1
-
c
1
b
2
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,
(5)
J
E
X
P
-
4
X
,
c
1
,
c
2
=
1
2
V
c
1
c
1
-
c
1
b
2
+
1
2
V
c
2
c
2
-
c
2
b
2
+
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,where X
=
x
1
0
,
x
2
0
,
x
3
0
,
w
0
,
x
1
1
,
x
2
1
,
x
3
1
,
w
1
T represents the initial x
1,2
,
3 and w at two time levels whose background values are denoted as X
b. The superscript “T” denotes the transpose. The background error covariance matrix is simply set as the observational error covariance matrix R
=
d
i
a
g
{
4.0
,
4.0
,
4.0
,
0.04
,
4.0
,
4.0
,
4.0
,
0.04
}. t
0 is the initial time within an MTW. N represents the number of observations within the MTW. H is the linearized observation operator which is d
i
a
g
{
1
,
1
,
1
,
1
,
0
,
0
,
0
,
0
} in this study. Y
o
=
x
1
o
,
x
2
o
,
x
3
o
,
w
o
T is the observation vector in which x
1,2
,
3
o and w
o, respectively, represent the observations of x
1,2
,
3 and w. c
1
b and c
2
b indicate the background values of c
1 and c
2, respectively. The background error variances of c
1 and c
2, that is, V
c
1 and V
c
2, are set to 0.00552 and 0.0552.The gradient of the cost function with respect to the parameter can be obtained using the adjoint method [34]. Due to the nonlinearity of the Lorenz-63 model, many local minima of the cost function exist in the space of parameter values when the MTW is relatively long, which may cause the standard quasi-Newton method (e.g., L-BGFS) [35] fail to find the optimal solution. In this study, we employ a limited memory bundle method (LMBM) [16, 36, 37] to implement the minimization of 4D-Var. Different from other methods [17, 38] for solving the multiple minima issue, LMBM is a hybrid of the variable metric bundle methods [39] and the limited memory variable metric methods [40]. This method is a solver of large-scale nonsmooth global optimization which neither needs to solve the time-consuming direction-finding issue in the standard bundle method nor needs to increase the number of stored subgradients when the dimension of the problem grows. The necessary input for LMBM algorithm includes the number of control variables, the maximum bundle dimension, the upper limit for maximum number of stored corrections, the maximum number of stored corrections, the first guess of control variables, the cost function value, the tolerance for the first or second termination criterion, and the maximum number of iterations [41]. Trial and error tests suggest that the last two inputs should be set to 10−5 and 200.In addition to the observational interval, the sensitivity of the 4D-Var CPE with respect to the length of MTW is also studied. We examined five MTWs: 0.5, 1.0, 2.0, 5.0, and 10 TUs. For all four data assimilation experiments, 4D-Var optimizes the initial fields when observations are available; that is, the analysis interval is the same as the observational interval. The absolute error of the optimized initial field relative to the truth is used to assess the quality of the 4D-Var CPE.


### 2.6. Implementation of EnKF CPE

In this study, a derivative of deterministic EnKF which does not need to perturb observations, the ensemble adjustment Kalman filter (EAKF) [11], is employed to implement the EnKF CPE.For an EnKF algorithm, the analysis solution includes the following two parts. One is the updated ensemble mean formulated as(6)
X
-
a
=
X
-
b
+
K
Y
o
-
H
X
-
b
,where X
-
a and X
-
b are the analysis and background of ensemble mean (X
-) of state vector X and K is the Kalman-gain matrix sampled by a finite ensemble. The other is the analysis solution of ensemble perturbations, which depends on the version of EnKF.When observation errors are assumed to be uncorrelated, EAKF can sequentially assimilate observations with the following two steps. First, the observational increment in one medium is computed as follows:(7)
Δ
y
i
=
r
2
r
2
+
σ
y
p
2
-
1
y
i
p
-
y
-
p
+
σ
y
p
2
r
2
+
σ
y
p
2
y
-
y
-
p
,where y
-
p represents the prior ensemble mean (i.e., the model estimate) of observation y; r and σ
y
p denote the standard deviation of observational errors and the prior standard deviation of y. The ith prior ensemble of y, y
i
p, is usually obtained through applying a linear interpolation to the prior ensemble of state variable. In this study, y represents one of x
1
o
,
x
2
o
,
x
3
o
,
w
o.Second step projects the observational increment onto related model variables and parameters. Realizations of this step for four EnKF CPE experiments are different. Similar to (2)–(5), here, we also give the linear regression formulas for four experiments as follows: (8)
E
X
P
-
1
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
(9)
E
X
P
-
2
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
2
,
i
=
c
o
v
c
2
,
y
σ
y
p
2
Δ
y
i
,
(10)
E
X
P
-
3
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
1
,
i
=
c
o
v
c
1
,
y
σ
y
p
2
Δ
y
i
,
(11)
E
X
P
-
4
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
1
,
i
=
c
o
v
c
1
,
y
σ
y
p
2
Δ
y
i
,
E
X
P
-
4
:
Δ
c
2
,
i
=
c
o
v
c
2
,
y
σ
y
p
2
Δ
y
i
,where Δ
x
i is the contribution of y to the model variable x for the ith ensemble member. In this study, x represents one of x
1
0
,
x
2
0
,
x
3
0
,
w
0
,
x
1
1
,
x
2
1
,
x
3
1
,
w
1 with superscripts “0” and “1” indicating time t
-
1 and time t. c
o
v
(
x
,
y
) denotes the error covariance between the prior ensemble of x and the model-estimated ensemble of y. It is worth mentioning that all experiments adopt the same observational increments and (8)–(11) are the main analysis equations of four experiments.Coupling parameter estimation in EnKF is a process similar to multivariate adjustment in state estimation for a nonobservable variable. Without dynamical support, the model parameter ensemble is easier to suffer the filter divergence than the model state. Thus, variance inflation is necessary for parameter estimations. The parameter variance inflation scheme adopted in this study is the same as that in the previous studies [23, 24, 26, 27, 42]; that is, when the current spread of the parameter ensemble is less than one-tenth of the initial spread, it will be inflated to this amount. For the state inflation scheme, although many sophisticated adaptively multiplicative inflation schemes are available [43–46], due to the low dimension of the simple coupled model and the perfect OSSE configuration here, we apply the static-multiplicative inflation scheme [47] that inflates the perturbations of model states with a constant factor. Since the coupled system has two characteristic time scales, two independent inflation factors are, respectively, applied to the chaotic atmospheric variables x
1,2
,
3 and the slab oceanic variable w. Through repeating the EnKF CPE experiment in the 2-dimensional space (
[
1
,
1.2
]
×
[
1
,
1.2
]
) of the two inflation factors with the interval of 0.01, the inflation factors that minimize the time-averaged (last 50 TUs) absolute errors of ensemble means of model states are chosen as the best pair to perform the state inflation.The ensemble initial conditions of model states are generated through adding a Gaussian noise to the initial fieldΠ. The initial ensembles of c
2 and c
1 are generated by adding a Gaussian noise with standard deviations of 0.055 and 0.0055 to the biased values 1.1 and 0.11. All four data assimilation experiments in the EnKF CPE start from the same ensemble initial conditions. The observations used in the EnKF CPE are the same as those used in the 4D-Var CPE. The absolute error of ensemble mean relative to the truth is used to evaluate the quality of the EnKF CPE. To investigate the dependence of the EnKF CPE on ensemble size, four ensemble sizes (10, 20, 50, and 100) are examined.


## 2.1. The Simple Coupled Model

To study coupling parameter estimation (CPE) implemented by 4D-Var and EnKF so as to answer the question (i.e., can 4D-Var or EnKF CPE improve representation of the model air-sea coupling process?) posed in Section1, we employ a simple “climate” model [25] which consists of a chaotic atmosphere model [29] coupled to a slab ocean model to simulate the interaction of the fast atmosphere and slow ocean:(1)
x
˙
1
=
-
σ
x
1
+
σ
x
2
,
x
˙
2
=
-
x
1
x
3
+
1
+
c
1
w
κ
x
1
-
x
2
,
x
˙
3
=
x
1
x
2
-
b
x
3
,
O
m
w
˙
=
c
2
x
2
-
O
d
w
+
S
m
+
S
s
c
o
s
2
π
t
S
p
d
,where an overdot denotes time tendency; x
1, x
2, and x
3 (hereafter denoted by x
1,2
,
3 if they are presented together) are the high-frequency variables of the atmosphere and w represents the slab ocean. (O
m, O
d) define the oceanic time scale, where O
m is the ratio of heat capacity between the slab ocean and the chaotic atmosphere while O
d is the damping coefficient of the temperature w (°C) of the slab ocean. O
m is a dimensionless parameter. The unit of O
d is 1
/
Δ
t, where Δ
t represents the dimensionless time step. The “atmosphere” and the “ocean” interact with each other through the coupling coefficients c
1 and c
2. Here, the unit of c
1 is (
C
°
)
-
1 while c
2 is a dimensionless parameter. The external solar forcing is represented by S
m
+
S
s
c
o
s
⁡
(
2
π
t
/
S
p
d
), where S
m represents the annual mean temperature (°C) of the slab ocean and S
s (in °C) and S
p
d (in Δ
t) represent the amplitude and period of the model seasonal cycle.While the detailed description of the model construction can be found in [25], here, we only comment on the setting of model parameter values. The “atmospheric” parameters σ, κ, and b take their standard values of 9.95, 29, and 8/3. For maintaining the chaotic nature of the “atmosphere” and the stability of the system, the values of c
1 and c
2 are chosen as 0.1 (
C
°
)
-
1 and 1. O
m and O
d are set to 10 and 1, which defines the oceanic time scale as 10 times of the atmospheric time scale. The parameters S
m and S
s are set as 10°C and 1°C, respectively. The S
p
d is chosen as 10 Time Units (TU, 1 TU = 100 time steps = 100Δ
t) so that the period of the forcing is comparable with the oceanic time scale, defining the time scale of the model seasonal cycle. Given the value of S
p
d, the model calendar year is defined as 10 TUs. If we assume that one year has 360 model days, Δ
t is equivalent to 0.36 model days.Using a leap-frog time stepping scheme with a Robert-Asselin time filter [31, 32] and starting from the initial conditions (
x
1
,
x
2
,
x
3
,
w
)
=
(
0
,
1
,
0
,
0
) (note that since the leap-frog time stepping scheme is applied, initial conditions at both time 0 and time 1 are set to (
0,1
,
0,0
).), the model is first freely run for 104 TUs by setting the parameters (σ, κ, b, c
1, c
2, O
m, O
d, S
m, S
s, S
p
d, and γ), where γ is the time filtering coefficient, to be the default values as prescribed in Table 1.Table 1
Default values and values used in OSSE of the model parameters.
Parameter
Default value
Value used in OSSE
Unit
σ
9.95
9.95
—
κ
29
29
—
b
8/3
8/3
—
c
1
0.1
0.11
(°C)−1
c
2
1
1.1
—
O
m
10
10
—
O
d
1
1
1/Δ
t
S
m
10
10
°C
S
s
1
1
°C
S
p
d
10
10
Δ
t
γ
0.25
0.25
—


## 2.2. OSSE Setup

OSSE that takes a Nature Run as the “truth” is an effective way to preliminarily study the validity of data assimilation algorithm and/or assess the observing system impact [30]. In this study, we setup an OSSE to study the EnKF CPE and the 4D-Var CPE. For simplicity, we assume that the model error only comes from the errors of the parameters to be estimated (see Table 1). The coupled model with the default parameter values produces the true solution for the parameter estimation problem. The “observations” are the samples of the “truth” model states after the spin-up of 104 TUs. A random Gaussian noise is superimposed on the “truth” values every 5 time steps for x
1,2
,
3 and every 20 time steps for w, respectively, to produce the corresponding observations. Note that, to simulate the real world sampling frequency, the atmospheric observations are taken more frequently (4 times) than the oceanic observations. The standard deviations of observational errors are 2 for x
1,2
,
3 and 0.2 for w, respectively.The parameters to be estimated arec
1 and c
2 that are the coupling parameters between the ocean and the atmosphere. Through examining the estimation of c
1 and c
2, we can directly study the impact of coupling parameter estimation on the representation of coupling processes. First guess of the parameters c
1 and c
2 is 0.11 (the “truth” is 0.1) and 1.1 (the “truth” is 1.0), respectively. From the initial condition (
x
1
,
x
2
,
x
3
,
w
)
=
(
0
,
1
,
0
,
0
), the assimilation model with the first guess of c
1 and c
2 is spun up for 104 TUs to produce the assimilation initial conditions of x
1,2
,
3 and w (denote as Π). Continuous integrations from these initial values of the assimilation model serve as a control run without any observational constraint. Therefore, this is a free model run, that is, without data assimilation.To study the 4D-Var CPE and EnKF CPE, we design four experiments for each algorithm:(
1
) EXP-1: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean while keeping c
1 and c
2 constant; (
2
) EXP-2: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean as well as c
2 while keeping c
1 constant; (
3
) EXP-3: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean as well as c
1 while keeping c
2 constant; (
4
) EXP-4: the atmospheric and oceanic observations are used to estimate the states of the atmosphere and the ocean, c
1 and c
2. Table 2 gives the details of assimilation schemes of these four experiments. The data assimilation period is set to be 100 TUs. Parameter estimation is activated after 10 TUs of state estimation which reaches quasi-equilibrium so that the state-parameter covariance used for parameter estimation is signal-dominated [25]. In addition, the leap-frog time stepping requires two time-level adjustments [33] (i.e., observations at time t are used to adjust the model states at time t and time t
-
1). And also, to investigate the sensitivity of the 4D-Var CPE and the EnKF CPE with respect to the observational interval, four pairs of (atmosphere and ocean) observational intervals: (0.05, 0.2), (0.1, 0.4), (0.2, 0.8), and (0.5, 2.0) TUs are examined. Finally, the results of the last 50 TUs are used for CPE evaluation. Note that, given the low-order model and the scalar nature of parameters, the absolute error of the estimated parameter or the air-sea coupling processes which refer to the truth is the main criterion of the CPE evaluation.Table 2
The instantaneous state and parameter estimations in four data assimilation experiments.
Experiment
Observations
Adjusted variables and parameters
EXP-1
x
1,2
,
3
o, w
o
x
1,2,3, w
EXP-2
x
1,2
,
3
o, w
o
x
1,2,3, w, c
2
EXP-3
x
1,2
,
3
o, w
o
x
1,2,3, w, c
1
EXP-4
x
1,2
,
3
o, w
o
x
1,2,3, w, c
1, c
2We see from time series ofx
2 (Figure 1(a)) and w (Figure 1(b)) in the truth, model control, and “observations” that the erroneous parameters can lead the model states of the free run to depart far away from the truth and observations. According to the governing equation (Equation (1)) of the simple coupled model, the Lorenz-atmosphere is forced by the term c
1
w from the slab ocean while the slab ocean couples with the atmosphere via the term c
2
x
2. Thus, we examine the c
1
w and c
2
x
2 errors to study the representation of the air-sea coupling processes in this study.Figure 1
Time series of (a)x
2 for the first 10 TUs and (b) w for 100 TUs of the assimilation period derived from the truth (red), model control run (blue), and observations (black solid circle), respectively.
(a)
(b)


## 2.3. Scale Analysis of the Perturbed Terms in the Model Equations

Within the OSSE framework, here, we make a scale analysis of the size of the perturbed terms in the model equations. Since the model error is assumed to arise from the uncertainties ofc
2 and c
1, in the ocean, the impact of an erroneous parameter c
2 and an erroneous atmosphere x
2 is contained in the term c
2
x
2. The error is δ
c
2
x
2
=
x
2
δ
c
2
+
c
2
δ
x
2. The characteristic scales of δ
c
2 and δ
x
2 are 0.1 and 10, respectively. Thus, given the default value (i.e., 1.0) of c
2 and the characteristic scale (i.e., 10) of x
2, the scale [i.e., the term x
2
δ
c
2 which is about 1.0] of correction of the oceanic parameter c
2 would be one order smaller than the scale [i.e., the term c
2
δ
x
2 which is about 10] of correction of the atmospheric state x
2.Similarly, in the atmosphere, the impact of an erroneous parameterc
1 and an erroneous ocean w is contained in the term c
1
w. The error is δ
c
1
w
=
w
δ
c
1
+
c
1
δ
w. The characteristic scales of δ
c
1 and δ
w are 0.01 and 1.0, respectively. Thus, given the default value (i.e., 0.1) of c
1 and the characteristic scale (i.e., 10) of w, the scale [i.e., the term w
δ
c
1 which is about 0.1] of correction of the atmospheric parameter c
1 would be the same as the scale [i.e., the term c
1
δ
w which is about 0.1] of correction of the oceanic state w.To summarize, the error of the oceanic state is one order smaller than the error of the atmospheric state. Thus, estimating the oceanic parameterc
2 will be more difficult than estimating the atmospheric parameter c
1.


## 2.4. Model Sensitivities with respect to Parameters

We perform a sensitivity study for seven empirical parameters (σ, κ, b, O
m, O
d, c
1, and c
2) using the same strategy as [26]. Each examined parameter is perturbed through adding a Gaussian noise to its default value (a percentage of the default value as the standard deviation of the noise) while other parameters remain unperturbed. With an ensemble size of 20, the model is freely run for 20 TUs. It is found that, as long as the perturbed parameter values can maintain the stochastic nature of the “atmosphere,” the examined ensemble spread does not have much dependence on the perturbation amplitude (i.e., the standard deviation of the above Gaussian noise) but the spin-up time is a little longer for a small perturbation. For a 5% perturbation (i.e., the standard deviation of the Gaussian noise is 5% of the default value), Figure 2 shows the time series of the ensemble spreads of x
2 (Figure 2(a)) and w (Figure 2(b)), where the black, blue, red, green, yellow, dashed-black, and pink curves, respectively, represent the results of b, c
1, c
2, κ, O
d, O
m, and σ. Here, the ensemble spreads of x
2 and w have been normalized by their own climatological standard deviations (14.5 for x
2 and 1.5 for w). We can see that different variables have different response times (defined as the spin-up period of ensemble spread) with regard to various parameters. For the four atmospheric parameters, the response times of x
2 and w are about 3 TUs and 6 TUs, respectively. For the three oceanic parameters, the response times of x
2 and w are about 5 TUs and 10 TUs, respectively. Therefore, if 4D-Var is used to perform CPE of the oceanic parameters, the typical lengths of MTWs for x
2 and w could, respectively, be 5 TUs and 10 TUs. For example, for x
2, an MTW less than 5 TUs is regarded as a short window while an MTW greater than 5 TUs is regarded as a long window.Figure 2
Time series of sensitivities of (a)x
2 and (b) w with respect to σ (pink), κ (green), b (black), O
m (dashed-black), O
d (yellow), c
1 (blue), and c
2 (red). Note that here the sensitivities of x
2 and w are computed as the ensemble spreads normalized by their climatological standard deviations (i.e., 14.5 and 1.5). Parameters σ, κ, b, O
m, and c
2 are dimensionless while the units of O
d and c
1 are 1
/
Δ
t and (°C)−1, where Δ
t is the dimensionless time step.
(a)
(b)From Figure2, we can see that, among four atmospheric parameters (σ, κ, b and c
1), both x
2 and w have the largest sensitivities with respect to κ (green) and c
1 (blue) during the spin-up period. Similarly, for the three oceanic parameters (O
d, O
m, and c
2), x
2 and w are most sensitive to the coupling parameter c
2 (red). Combining the question presented in Section 1, two coupling parameters c
1 and c
2 are chosen as the parameters to be estimated in the following CPE experiments.


## 2.5. Implementation of 4D-Var CPE

Consistent with Table2, the formulas of the cost functions for four experiments are(2)
J
E
X
P
-
1
X
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
R
-
1
H
X
t
i
-
Y
o
t
i
,
(3)
J
E
X
P
-
2
X
,
c
2
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
V
c
2
c
2
-
c
2
b
2
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,
(4)
J
E
X
P
-
3
X
,
c
1
=
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
V
c
1
c
1
-
c
1
b
2
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,
(5)
J
E
X
P
-
4
X
,
c
1
,
c
2
=
1
2
V
c
1
c
1
-
c
1
b
2
+
1
2
V
c
2
c
2
-
c
2
b
2
+
1
2
X
t
0
-
X
b
T
B
-
1
X
t
0
-
X
b
+
1
2
∑
i
=
1
N
H
X
t
i
-
Y
o
t
i
T
·
R
-
1
H
X
t
i
-
Y
o
t
i
,where X
=
x
1
0
,
x
2
0
,
x
3
0
,
w
0
,
x
1
1
,
x
2
1
,
x
3
1
,
w
1
T represents the initial x
1,2
,
3 and w at two time levels whose background values are denoted as X
b. The superscript “T” denotes the transpose. The background error covariance matrix is simply set as the observational error covariance matrix R
=
d
i
a
g
{
4.0
,
4.0
,
4.0
,
0.04
,
4.0
,
4.0
,
4.0
,
0.04
}. t
0 is the initial time within an MTW. N represents the number of observations within the MTW. H is the linearized observation operator which is d
i
a
g
{
1
,
1
,
1
,
1
,
0
,
0
,
0
,
0
} in this study. Y
o
=
x
1
o
,
x
2
o
,
x
3
o
,
w
o
T is the observation vector in which x
1,2
,
3
o and w
o, respectively, represent the observations of x
1,2
,
3 and w. c
1
b and c
2
b indicate the background values of c
1 and c
2, respectively. The background error variances of c
1 and c
2, that is, V
c
1 and V
c
2, are set to 0.00552 and 0.0552.The gradient of the cost function with respect to the parameter can be obtained using the adjoint method [34]. Due to the nonlinearity of the Lorenz-63 model, many local minima of the cost function exist in the space of parameter values when the MTW is relatively long, which may cause the standard quasi-Newton method (e.g., L-BGFS) [35] fail to find the optimal solution. In this study, we employ a limited memory bundle method (LMBM) [16, 36, 37] to implement the minimization of 4D-Var. Different from other methods [17, 38] for solving the multiple minima issue, LMBM is a hybrid of the variable metric bundle methods [39] and the limited memory variable metric methods [40]. This method is a solver of large-scale nonsmooth global optimization which neither needs to solve the time-consuming direction-finding issue in the standard bundle method nor needs to increase the number of stored subgradients when the dimension of the problem grows. The necessary input for LMBM algorithm includes the number of control variables, the maximum bundle dimension, the upper limit for maximum number of stored corrections, the maximum number of stored corrections, the first guess of control variables, the cost function value, the tolerance for the first or second termination criterion, and the maximum number of iterations [41]. Trial and error tests suggest that the last two inputs should be set to 10−5 and 200.In addition to the observational interval, the sensitivity of the 4D-Var CPE with respect to the length of MTW is also studied. We examined five MTWs: 0.5, 1.0, 2.0, 5.0, and 10 TUs. For all four data assimilation experiments, 4D-Var optimizes the initial fields when observations are available; that is, the analysis interval is the same as the observational interval. The absolute error of the optimized initial field relative to the truth is used to assess the quality of the 4D-Var CPE.


## 2.6. Implementation of EnKF CPE

In this study, a derivative of deterministic EnKF which does not need to perturb observations, the ensemble adjustment Kalman filter (EAKF) [11], is employed to implement the EnKF CPE.For an EnKF algorithm, the analysis solution includes the following two parts. One is the updated ensemble mean formulated as(6)
X
-
a
=
X
-
b
+
K
Y
o
-
H
X
-
b
,where X
-
a and X
-
b are the analysis and background of ensemble mean (X
-) of state vector X and K is the Kalman-gain matrix sampled by a finite ensemble. The other is the analysis solution of ensemble perturbations, which depends on the version of EnKF.When observation errors are assumed to be uncorrelated, EAKF can sequentially assimilate observations with the following two steps. First, the observational increment in one medium is computed as follows:(7)
Δ
y
i
=
r
2
r
2
+
σ
y
p
2
-
1
y
i
p
-
y
-
p
+
σ
y
p
2
r
2
+
σ
y
p
2
y
-
y
-
p
,where y
-
p represents the prior ensemble mean (i.e., the model estimate) of observation y; r and σ
y
p denote the standard deviation of observational errors and the prior standard deviation of y. The ith prior ensemble of y, y
i
p, is usually obtained through applying a linear interpolation to the prior ensemble of state variable. In this study, y represents one of x
1
o
,
x
2
o
,
x
3
o
,
w
o.Second step projects the observational increment onto related model variables and parameters. Realizations of this step for four EnKF CPE experiments are different. Similar to (2)–(5), here, we also give the linear regression formulas for four experiments as follows: (8)
E
X
P
-
1
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
(9)
E
X
P
-
2
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
2
,
i
=
c
o
v
c
2
,
y
σ
y
p
2
Δ
y
i
,
(10)
E
X
P
-
3
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
1
,
i
=
c
o
v
c
1
,
y
σ
y
p
2
Δ
y
i
,
(11)
E
X
P
-
4
:
Δ
x
i
=
c
o
v
x
,
y
σ
y
p
2
Δ
y
i
,
Δ
c
1
,
i
=
c
o
v
c
1
,
y
σ
y
p
2
Δ
y
i
,
E
X
P
-
4
:
Δ
c
2
,
i
=
c
o
v
c
2
,
y
σ
y
p
2
Δ
y
i
,where Δ
x
i is the contribution of y to the model variable x for the ith ensemble member. In this study, x represents one of x
1
0
,
x
2
0
,
x
3
0
,
w
0
,
x
1
1
,
x
2
1
,
x
3
1
,
w
1 with superscripts “0” and “1” indicating time t
-
1 and time t. c
o
v
(
x
,
y
) denotes the error covariance between the prior ensemble of x and the model-estimated ensemble of y. It is worth mentioning that all experiments adopt the same observational increments and (8)–(11) are the main analysis equations of four experiments.Coupling parameter estimation in EnKF is a process similar to multivariate adjustment in state estimation for a nonobservable variable. Without dynamical support, the model parameter ensemble is easier to suffer the filter divergence than the model state. Thus, variance inflation is necessary for parameter estimations. The parameter variance inflation scheme adopted in this study is the same as that in the previous studies [23, 24, 26, 27, 42]; that is, when the current spread of the parameter ensemble is less than one-tenth of the initial spread, it will be inflated to this amount. For the state inflation scheme, although many sophisticated adaptively multiplicative inflation schemes are available [43–46], due to the low dimension of the simple coupled model and the perfect OSSE configuration here, we apply the static-multiplicative inflation scheme [47] that inflates the perturbations of model states with a constant factor. Since the coupled system has two characteristic time scales, two independent inflation factors are, respectively, applied to the chaotic atmospheric variables x
1,2
,
3 and the slab oceanic variable w. Through repeating the EnKF CPE experiment in the 2-dimensional space (
[
1
,
1.2
]
×
[
1
,
1.2
]
) of the two inflation factors with the interval of 0.01, the inflation factors that minimize the time-averaged (last 50 TUs) absolute errors of ensemble means of model states are chosen as the best pair to perform the state inflation.The ensemble initial conditions of model states are generated through adding a Gaussian noise to the initial fieldΠ. The initial ensembles of c
2 and c
1 are generated by adding a Gaussian noise with standard deviations of 0.055 and 0.0055 to the biased values 1.1 and 0.11. All four data assimilation experiments in the EnKF CPE start from the same ensemble initial conditions. The observations used in the EnKF CPE are the same as those used in the 4D-Var CPE. The absolute error of ensemble mean relative to the truth is used to evaluate the quality of the EnKF CPE. To investigate the dependence of the EnKF CPE on ensemble size, four ensemble sizes (10, 20, 50, and 100) are examined.


## 3. Results of 4D-Var CPE

In this section, we first evaluate the quality of the 4D-Var-estimated coupling parameters and then investigate the impact of the estimated-parameters on the representation of coupling processes.


### 3.1. Coupling Parameters

Figure3 shows the dependence of time mean normalized absolute errors of c
2 (Figures 3(a) and 3(b)) and c
1 (Figures 3(c) and 3(d)) on the observational interval and the length of MTW for EXP-2 (Figure 3(a)), EXP-3 (Figure 3(c)), and EXP-4 (Figures 3(b) and 3(d)). Note that the absolute errors of c
2 and c
1 here are normalized by their initial errors (0.1 and 0.01). Note that we use observational interval instead of atmospheric observational interval in the following text. In EXP-4, 4D-Var can well estimate c
1 for all observational intervals and MTWs. For short MTWs (like 50 and 100 time steps), 4D-Var cannot effectively estimate c
2 with too small (like 5 and 10 time steps) or too large (like 50 time steps) observational intervals. According to the analysis in Section 2.3, it is more difficult to estimate c
2 than c
1. For short MTWs, no sufficient signals of c
2 are implied in the cost function of 4D-Var, which may increase the possibility of the failure of parameter estimation. In EXP-2 that only estimates c
2, we can see that 4D-Var fails in most cases except the situations with short MTWs (like 50 and 100 time steps) and a moderate observational interval (like 20 time steps). Since the model error is mainly caused by the error of c
1, 4D-Var cannot accurately estimate c
2 without estimation of c
1. In contrast, 4D-Var can effectively estimate c
1 in most cases in EXP-3. Once c
1 is estimated, the signal-to-noise ratio of c
2 in the cost function is correspondingly enhanced, which accordingly improves the quality of the estimated c
2 (see Figures 3(a) and 3(b)). Additionally, although contribution of c
2 to the model error is much less than c
1, introducing the estimation of c
2 to EXP-3 still can more or less further improve the accuracy of the estimated c
1 (see Figures 3(c) and 3(d)).Figure 3
Time-averaged normalized absolute errors of (a)c
2 in EXP-2, (b) c
2 in EXP-4, (c) c
1 in EXP-3, and (d) c
1 in EXP-4 of the 4D-Var CPE. The absolute errors of c
2 and c
1 are normalized by their initial errors, that is, 0.1 and 0.01, respectively. The x-axis is the atmospheric observational interval while the y-axis denotes the minimization time window (MTW) in 4D-Var. The black solid circles stand for the scatter points used to create the contours. The dashed curve represents the 1.0 contour. Panels (a) and (b) use the upper shade scale while panels (c) and (d) use the lower shade scale. Note that 1 Time Unit = 100 time steps. EXP-1 only performs state estimation. EXP-2 instantaneously conducts state estimation and parameter estimation of c
2. EXP-3 instantaneously conducts state estimation and parameter estimation of c
1. EXP-4 instantaneously conducts state estimation and parameter estimations of c
1 and c
2.
(a)
(b)
(c)
(d)From Figures3(b) and 3(d), we can also see that the analysis accuracies of two parameters (especially for c
2) strongly depend on both MTW and observational intervals when MTW is less than 500 time steps. This means that, to obtain good analysis solutions of both c
1 and c
2, an appropriate MTW (about 500 time steps here) should be used. Note that, due to the limited number of experiments performed in this study, for example, no experiment with MTW between 200 and 500 time steps is conducted; here, the appropriate MTW is an approximation.


### 3.2. Coupling Processes

Figure4 displays the time-averaged normalized absolute errors of the sea-to-air coupling process (represented by c
1
w) for EXP-1 (Figure 4(a)), EXP-2 (Figure 4(b)), EXP-3 (Figure 4(c)), and EXP4 (Figure 4(d)). Here, the absolute error is normalized by the climatological standard deviation (0.15) of c
1
w. Compared to EXP-1, estimating c
2 obtains nearly the same analysis of c
1
w (see Figure 4(b)). However, once c
1 is estimated, the sea-to-air coupling process is significantly improved (compare Figure 4(a) with Figure 4(c)). The given parameter c
1 has been suitably corrected (see Figure 3(c)) in EXP-3, further estimating that c
2 can also partly improve the sea-to-air coupling process (compare Figure 4(c) with Figure 4(d)).Figure 4
Time-averaged normalized absolute errors ofc
1
w in (a) EXP-1, (b) EXP-2, (c) EXP-3, and (d) EXP-4 of the 4D-Var CPE. The absolute error of c
1
w is normalized by its climatological standard deviation (0.15). All panels use the same shade scale.
(a)
(b)
(c)
(d)Figure5 shows the same results as Figure 4 but for the air-to-sea coupling process (represented by c
2
x
2). Based on EXP-1, estimating c
2 cannot significantly improve the quality of c
2
x
2 (compare Figure 5(b) with Figure 5(a)). However, when c
1 is further estimated, the quality of the air-to-sea coupling processes is greatly enhanced, especially for long MTWs. In addition, comparison between Figure 5(c) and Figure 5(d) justifies that given c
1 has been corrected and the quality of the air-to-sea coupling processes can also be somehow enhanced by estimating c
2. Thus, unlike the sea-to-air coupling process (c
1
w) which is governed by the c
1 accuracy, the air-to-sea coupling process (c
2
x
2) is affected by both c
1 and c
2 greatly.Figure 5
Same as Figure4 but for c
2
x
2. The absolute error of c
2
x
2 is normalized by its climatological standard deviation (14.5).
(a)
(b)
(c)
(d)To look into the detailed performance of 4D-Var CPE, we choose an appropriate MTW (i.e., 500 time steps) and a moderate observational interval (i.e., 10 time steps) to examine the time series of the absolute errors of the air-to-sea (Figure6(b)) and the sea-to-air (Figure 6(a)) coupling processes for EXP-1 (red) and EXP-4 (blue). Relative to the state estimation, coupling parameter estimation of c
1 and c
2 can markedly enhance the accuracy of coupling processes. According to Section 2.4, model states are more sensitive to c
1 than c
2. Scale analysis (Section 2.3) also justifies that estimating c
2 is more difficult than estimating c
1. Thus, the model error in this study is mainly attributed to the bias of c
1, causing that the improvement of the air-to-sea coupling process is less than that of the sea-to-air coupling process.Figure 6
Time series (20–90 TUs) of the absolute errors of (a)c
1
w and (b) c
2
x
2 in the 4D-Var CPE with the observational interval and the minimization time window being 10 and 500 time steps. Here, the red and blue curves represent the results in EXP-1 and EXP-4, respectively.
(a)
(b)Analyses above show that, compared to traditional state estimation, 4D-Var CPE can improve the representation of air-sea coupling processes through coupling parameter estimation.


## 3.1. Coupling Parameters

Figure3 shows the dependence of time mean normalized absolute errors of c
2 (Figures 3(a) and 3(b)) and c
1 (Figures 3(c) and 3(d)) on the observational interval and the length of MTW for EXP-2 (Figure 3(a)), EXP-3 (Figure 3(c)), and EXP-4 (Figures 3(b) and 3(d)). Note that the absolute errors of c
2 and c
1 here are normalized by their initial errors (0.1 and 0.01). Note that we use observational interval instead of atmospheric observational interval in the following text. In EXP-4, 4D-Var can well estimate c
1 for all observational intervals and MTWs. For short MTWs (like 50 and 100 time steps), 4D-Var cannot effectively estimate c
2 with too small (like 5 and 10 time steps) or too large (like 50 time steps) observational intervals. According to the analysis in Section 2.3, it is more difficult to estimate c
2 than c
1. For short MTWs, no sufficient signals of c
2 are implied in the cost function of 4D-Var, which may increase the possibility of the failure of parameter estimation. In EXP-2 that only estimates c
2, we can see that 4D-Var fails in most cases except the situations with short MTWs (like 50 and 100 time steps) and a moderate observational interval (like 20 time steps). Since the model error is mainly caused by the error of c
1, 4D-Var cannot accurately estimate c
2 without estimation of c
1. In contrast, 4D-Var can effectively estimate c
1 in most cases in EXP-3. Once c
1 is estimated, the signal-to-noise ratio of c
2 in the cost function is correspondingly enhanced, which accordingly improves the quality of the estimated c
2 (see Figures 3(a) and 3(b)). Additionally, although contribution of c
2 to the model error is much less than c
1, introducing the estimation of c
2 to EXP-3 still can more or less further improve the accuracy of the estimated c
1 (see Figures 3(c) and 3(d)).Figure 3
Time-averaged normalized absolute errors of (a)c
2 in EXP-2, (b) c
2 in EXP-4, (c) c
1 in EXP-3, and (d) c
1 in EXP-4 of the 4D-Var CPE. The absolute errors of c
2 and c
1 are normalized by their initial errors, that is, 0.1 and 0.01, respectively. The x-axis is the atmospheric observational interval while the y-axis denotes the minimization time window (MTW) in 4D-Var. The black solid circles stand for the scatter points used to create the contours. The dashed curve represents the 1.0 contour. Panels (a) and (b) use the upper shade scale while panels (c) and (d) use the lower shade scale. Note that 1 Time Unit = 100 time steps. EXP-1 only performs state estimation. EXP-2 instantaneously conducts state estimation and parameter estimation of c
2. EXP-3 instantaneously conducts state estimation and parameter estimation of c
1. EXP-4 instantaneously conducts state estimation and parameter estimations of c
1 and c
2.
(a)
(b)
(c)
(d)From Figures3(b) and 3(d), we can also see that the analysis accuracies of two parameters (especially for c
2) strongly depend on both MTW and observational intervals when MTW is less than 500 time steps. This means that, to obtain good analysis solutions of both c
1 and c
2, an appropriate MTW (about 500 time steps here) should be used. Note that, due to the limited number of experiments performed in this study, for example, no experiment with MTW between 200 and 500 time steps is conducted; here, the appropriate MTW is an approximation.


## 3.2. Coupling Processes

Figure4 displays the time-averaged normalized absolute errors of the sea-to-air coupling process (represented by c
1
w) for EXP-1 (Figure 4(a)), EXP-2 (Figure 4(b)), EXP-3 (Figure 4(c)), and EXP4 (Figure 4(d)). Here, the absolute error is normalized by the climatological standard deviation (0.15) of c
1
w. Compared to EXP-1, estimating c
2 obtains nearly the same analysis of c
1
w (see Figure 4(b)). However, once c
1 is estimated, the sea-to-air coupling process is significantly improved (compare Figure 4(a) with Figure 4(c)). The given parameter c
1 has been suitably corrected (see Figure 3(c)) in EXP-3, further estimating that c
2 can also partly improve the sea-to-air coupling process (compare Figure 4(c) with Figure 4(d)).Figure 4
Time-averaged normalized absolute errors ofc
1
w in (a) EXP-1, (b) EXP-2, (c) EXP-3, and (d) EXP-4 of the 4D-Var CPE. The absolute error of c
1
w is normalized by its climatological standard deviation (0.15). All panels use the same shade scale.
(a)
(b)
(c)
(d)Figure5 shows the same results as Figure 4 but for the air-to-sea coupling process (represented by c
2
x
2). Based on EXP-1, estimating c
2 cannot significantly improve the quality of c
2
x
2 (compare Figure 5(b) with Figure 5(a)). However, when c
1 is further estimated, the quality of the air-to-sea coupling processes is greatly enhanced, especially for long MTWs. In addition, comparison between Figure 5(c) and Figure 5(d) justifies that given c
1 has been corrected and the quality of the air-to-sea coupling processes can also be somehow enhanced by estimating c
2. Thus, unlike the sea-to-air coupling process (c
1
w) which is governed by the c
1 accuracy, the air-to-sea coupling process (c
2
x
2) is affected by both c
1 and c
2 greatly.Figure 5
Same as Figure4 but for c
2
x
2. The absolute error of c
2
x
2 is normalized by its climatological standard deviation (14.5).
(a)
(b)
(c)
(d)To look into the detailed performance of 4D-Var CPE, we choose an appropriate MTW (i.e., 500 time steps) and a moderate observational interval (i.e., 10 time steps) to examine the time series of the absolute errors of the air-to-sea (Figure6(b)) and the sea-to-air (Figure 6(a)) coupling processes for EXP-1 (red) and EXP-4 (blue). Relative to the state estimation, coupling parameter estimation of c
1 and c
2 can markedly enhance the accuracy of coupling processes. According to Section 2.4, model states are more sensitive to c
1 than c
2. Scale analysis (Section 2.3) also justifies that estimating c
2 is more difficult than estimating c
1. Thus, the model error in this study is mainly attributed to the bias of c
1, causing that the improvement of the air-to-sea coupling process is less than that of the sea-to-air coupling process.Figure 6
Time series (20–90 TUs) of the absolute errors of (a)c
1
w and (b) c
2
x
2 in the 4D-Var CPE with the observational interval and the minimization time window being 10 and 500 time steps. Here, the red and blue curves represent the results in EXP-1 and EXP-4, respectively.
(a)
(b)Analyses above show that, compared to traditional state estimation, 4D-Var CPE can improve the representation of air-sea coupling processes through coupling parameter estimation.


## 4.  Results of EnKF CPE

In this section, same as in the 4D-Var case, we first assess the quality of EnKF-estimated coupling parameters and then investigate the impact of the estimated-parameters on the representation of coupling processes.


### 4.1. Coupling Parameters

Figure7 gives the same results as Figure 3 but for EnKF CPE. Here, the y-axis represents ensemble size. Different from 4D-Var (Figure 3(a)), EnKF can effectively estimate c
2 in EXP-2 (Figure 7(a)). Based on EXP-2, EXP-4 can further reduce the error of estimated c
2 through adjusting c
1. Both EXP-3 and EXP-4 can do a good job for the estimation of c
1. Thus, it is not so necessary that EnKF should firstly estimate c
1 rather than c
2. This may be attributed to the high signal-to-noise ratios of c
2 implied in the error covariance between observations and c
2.Figure 7
Same as Figure3 but for the EnKF CPE. Note that panels (a) and (b) use the upper shade scale while panels (c) and (d) use the lower shade scale.
(a)
(b)
(c)
(d)


### 4.2. Coupling Processes

Figure8 plots the same results as Figure 4 but for EnKF CPE. From Figures 8(a) and 8(b), we can see that estimation of c
2 cannot significantly improve the representation of the sea-to-air coupling process. This is because the term c
1
w is directly controlled by c
1 and w rather than c
2. If we compare Figure 8(c) to Figure 8(a), we can see that estimating c
1 can substantially enhance the accuracy of c
1
w. When c
2 is further estimated based on EXP-3, the error of c
1
w is also further reduced.Figure 8
Same as Figure4 but for the EnKF CPE.
(a)
(b)
(c)
(d)Figure9 shows the same results as Figure 8 but for the air-to-sea coupling process. Relative to EXP-1, both EXP-2 and EXP-3 can reduce the error of c
2
x
2 while simultaneously estimating c
1 and c
2 produces the best results. Additionally, due to the low-order model in this study, the dependence of the quality of estimated air-to-sea coupling process on observational interval is stronger than that on ensemble size.Figure 9
Same as Figure5 but for the EnKF CPE.
(a)
(b)
(c)
(d)To look into the detailed performance of EnKF CPE, we set observational interval as 10 time steps and ensemble size as 20 to examine the time series of absolute errors of the air-to-sea (Figure10(b)) and the sea-to-air (Figure 10(a)) coupling processes for EXP-1 (red) and EXP-4 (blue). Compared to the air-to-sea coupling process, the improvement of the sea-to-air coupling process from EXP-1 to EXP-4 is more significant (i.e., the error of c
1
w is reduced by 92%, from 0.53 to 0.04).Figure 10
Same as Figure6 but for the EnKF CPE with 20-ensemble size.
(a)
(b)Analyses above show that EnKF CPE can also enhance the representation of the air-sea coupling processes through coupling parameter estimation relative to traditional state estimation.


## 4.1. Coupling Parameters

Figure7 gives the same results as Figure 3 but for EnKF CPE. Here, the y-axis represents ensemble size. Different from 4D-Var (Figure 3(a)), EnKF can effectively estimate c
2 in EXP-2 (Figure 7(a)). Based on EXP-2, EXP-4 can further reduce the error of estimated c
2 through adjusting c
1. Both EXP-3 and EXP-4 can do a good job for the estimation of c
1. Thus, it is not so necessary that EnKF should firstly estimate c
1 rather than c
2. This may be attributed to the high signal-to-noise ratios of c
2 implied in the error covariance between observations and c
2.Figure 7
Same as Figure3 but for the EnKF CPE. Note that panels (a) and (b) use the upper shade scale while panels (c) and (d) use the lower shade scale.
(a)
(b)
(c)
(d)


## 4.2. Coupling Processes

Figure8 plots the same results as Figure 4 but for EnKF CPE. From Figures 8(a) and 8(b), we can see that estimation of c
2 cannot significantly improve the representation of the sea-to-air coupling process. This is because the term c
1
w is directly controlled by c
1 and w rather than c
2. If we compare Figure 8(c) to Figure 8(a), we can see that estimating c
1 can substantially enhance the accuracy of c
1
w. When c
2 is further estimated based on EXP-3, the error of c
1
w is also further reduced.Figure 8
Same as Figure4 but for the EnKF CPE.
(a)
(b)
(c)
(d)Figure9 shows the same results as Figure 8 but for the air-to-sea coupling process. Relative to EXP-1, both EXP-2 and EXP-3 can reduce the error of c
2
x
2 while simultaneously estimating c
1 and c
2 produces the best results. Additionally, due to the low-order model in this study, the dependence of the quality of estimated air-to-sea coupling process on observational interval is stronger than that on ensemble size.Figure 9
Same as Figure5 but for the EnKF CPE.
(a)
(b)
(c)
(d)To look into the detailed performance of EnKF CPE, we set observational interval as 10 time steps and ensemble size as 20 to examine the time series of absolute errors of the air-to-sea (Figure10(b)) and the sea-to-air (Figure 10(a)) coupling processes for EXP-1 (red) and EXP-4 (blue). Compared to the air-to-sea coupling process, the improvement of the sea-to-air coupling process from EXP-1 to EXP-4 is more significant (i.e., the error of c
1
w is reduced by 92%, from 0.53 to 0.04).Figure 10
Same as Figure6 but for the EnKF CPE with 20-ensemble size.
(a)
(b)Analyses above show that EnKF CPE can also enhance the representation of the air-sea coupling processes through coupling parameter estimation relative to traditional state estimation.


## 5. Different Performances of EnKF CPE and 4D-Var CPE

Although both 4D-Var and EnKF are derived from Bayes’ rule, their implementations are different. In practice, 4D-Var attempts to obtain an optimal initial condition using the observational information within a specific MTW for each analysis step. The length of MTW is usually longer than the analysis interval, leading the observations being used more than once. In contrast, EnKF does not repeatedly use observations in this study (although it could do so in practice). From this point, it seems unfair to compare the performances of two algorithms. However, additional experiments for 4D-Var which only uses each observation for one time obtain worse results. Thus, we still can roughly compare the performances of EnKF CPE and 4D-Var CPE in this section.


### 5.1. Coupling Parameters

Comparison between Figure3(a) and Figure 7(a) demonstrates that EnKF-estimated c
2 is more accurate than that estimated by 4D-Var. Without estimating c
1, the signal implied in the error covariance between observation and c
2 is stronger than that implied in the cost function of 4D-Var. Since the model error is mainly attributed to the bias of c
1, the signal in the cost function of 4D-Var is dominated by c
1 rather than c
2. Under this circumstance, it is hard for 4D-Var to correctly estimate c
2. In EXP-3, results of EnKF CPE are still better than those for 4D-Var CPE which reduces the initial error of c
1 by 74% on average. In EXP-4, to facilitate the comparison, we use the same shade scale for Figure 3(b) and Figure 7(b) and Figure 3(d) and Figure 7(d). For both parameters, the performance of EnKF CPE with a moderate ensemble size (like 50) is comparable to that of 4D-Var CPE with an appropriate MTW (like 500 time steps).To summarize, EnKF is better than 4D-Var when a single parameter is estimated. The ensemble-evaluated error covariance is more effective than the cost function of 4D-Var.


### 5.2. Coupling Processes

We first attempt to compare the results of EXP-1 for 4D-Var CPE and EnKF CPE. For the sea-to-air coupling process represented byc
1
w, the quality of the estimated c
1
w is determined by the quality of w since c
1 is not estimated. From Figures 4(a) and 8(a), it seems that c
1
w produced by 4D-Var CPE is much better than that generated by EnKF CPE. Thus, we may speculate that the 4D-Var-estimated w is also better than EnKF-estimated w. However, results (not shown) demonstrate that the EnKF-produced w is much better than that produced by 4D-Var.To understand why a good EnKF-producedw leads a bad c
1
w, we examine an experiment for details with the observational interval, ensemble size, and MTW being 5 time steps, 50 time steps, and 500 time steps. Figure 11 plots the time series of w (Figure 11(a)) and c
1
w (Figure 11(b)) for 4D-Var CPE (black), truth (red), and EnKF CPE (blue) in EXP-1. Obviously, the EnKF-produced w is better than 4D-Var-produced w. We start from the following formulation of the absolute error of c
1
w to answer the question at the beginning of this paragraph:(12)
c
1
e
s
t
w
e
s
t
-
c
1
t
r
u
w
t
r
u
,where the superscripts “est” and “tru” indicate the estimation and truth of c
1 and w. As we stated in Section 2.2, the truth value of c
1 is 0.1 while the initial biased value of c
1 in 4D-Var CPE is 0.11 and ensemble means of c
1 for all ensemble sizes in EnKF CPE also approximate 0.11. Therefore, the computational process of (12) can be described as first scaling the truth value of w by 0.1 and the values of estimated w values by 0.11, then computing the absolute value of the difference. From Figure 11, it happens that the first step converts the worse 4D-Var-produced w to the appearing better c
1
w. Therefore, it is difficult to determine which algorithm produces better coupling processes without parameter estimation.Figure 11
Time series ofw (panel (a)) and c
1
w (panel (b)) for 4D-Var CPE (black) with MTW being 500 time steps, truth (red), and EnKF CPE (blue) with ensemble size being 50 in EXP-4. Here, the observational interval is set to 5 time steps.
(a)
(b)To fairly compare the performances of two algorithms, we focus on the results in EXP-4. From Figures4(d) and 8(d), we find that the representation of the sea-to-air coupling process produced by 4D-Var CPE with a MTW over 500 time steps is comparable to EnKF CPE. However, the EnKF CPE with ensemble size less than 50 improves the assimilation quality of 4D-Var CPE with a MTW shorter than 500 time steps by 59% (from 0.157 to 0.065) on average. Here, we choose an extreme case (i.e., the observational interval, ensemble size, and MTW being 50 time steps, 10 time steps, and 50 time steps) to compare the performance of two algorithms. Figure 12 shows time series of absolute errors of c
1
w for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. The EnKF CPE is better than the 4D-Var CPE. According to Figures 3(d) and 7(d), in the above parameter setting, the EnKF-estimated c
1 is much better than the 4D-Var-estimated c
1, since the MTW is too short to exactly retrieve c
1 for 4D-Var. In addition, we found that the 4D-Var-produced w is comparable to that produced by EnKF (not shown). Thus, the superiority of EnKF CPE over 4D-Var CPE is mainly attributed to the estimation of c
1.Figure 12
Time series of absolute errors ofc
1
w for 4D-Var CPE (blue) with MTW being 50 time steps and EnKF CPE (red) with 10-ensemble size in EXP-4. Here, the observational interval is 50 time steps.For the air-to-sea coupling process, comparison between Figure5(d) and Figure 9(d) reveals that 4D-Var CPE outperforms EnKF CPE when MTW exceeds a critical value (here is about 500 time steps) while EnKF CPE has advantages over 4D-Var CPE with a short MTW for a small observational interval. Here, we take two extreme cases to compare the performances of 4D-Var CPE and EnKF CPE. First, we check the results of two methods for the observational interval, ensemble size, and MTW being 50 time steps, 100 time steps, and 1000 time steps. Figure 13(a) plots time series of absolute errors of c
2
x
2 for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. After the spin-up period (about 40 TUs), the air-to-sea coupling process represented by 4D-Var CPE is better than that produced by EnKF CPE. Although the observational interval is very large and both 4D-Var CPE and EnKF CPE can effectively estimate c
2 (see Figures 3(b) and 7(b)), 4D-Var CPE can produce better x
2 than EnKF CPE (not shown). Second, we check the results (Figure 13(b)) of two methods for the observational interval, ensemble size, and MTW being 5 time steps, 10 time steps, and 50 time steps. According to Figures 3(b) and 7(b), it is difficult for 4D-Var CPE to estimate c
2 with a short MTW while EnKF CPE can well estimate c
2 with an ensemble size of 10. For x
2, 4D-Var CPE can produce similar results as EnKF CPE (not shown). Therefore, the EnKF CPE-produced c
2
x
2 is better than that produced by 4D-Var CPE.Figure 13
Time series of absolute errors ofc
2
x
2 for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. In panel (a) (panel (b)), values of the observational interval, ensemble size, and MTW are 50 (5) time steps, 100 (10) time steps, and 1000 (50) time steps.
(a)
(b)


## 5.1. Coupling Parameters

Comparison between Figure3(a) and Figure 7(a) demonstrates that EnKF-estimated c
2 is more accurate than that estimated by 4D-Var. Without estimating c
1, the signal implied in the error covariance between observation and c
2 is stronger than that implied in the cost function of 4D-Var. Since the model error is mainly attributed to the bias of c
1, the signal in the cost function of 4D-Var is dominated by c
1 rather than c
2. Under this circumstance, it is hard for 4D-Var to correctly estimate c
2. In EXP-3, results of EnKF CPE are still better than those for 4D-Var CPE which reduces the initial error of c
1 by 74% on average. In EXP-4, to facilitate the comparison, we use the same shade scale for Figure 3(b) and Figure 7(b) and Figure 3(d) and Figure 7(d). For both parameters, the performance of EnKF CPE with a moderate ensemble size (like 50) is comparable to that of 4D-Var CPE with an appropriate MTW (like 500 time steps).To summarize, EnKF is better than 4D-Var when a single parameter is estimated. The ensemble-evaluated error covariance is more effective than the cost function of 4D-Var.


## 5.2. Coupling Processes

We first attempt to compare the results of EXP-1 for 4D-Var CPE and EnKF CPE. For the sea-to-air coupling process represented byc
1
w, the quality of the estimated c
1
w is determined by the quality of w since c
1 is not estimated. From Figures 4(a) and 8(a), it seems that c
1
w produced by 4D-Var CPE is much better than that generated by EnKF CPE. Thus, we may speculate that the 4D-Var-estimated w is also better than EnKF-estimated w. However, results (not shown) demonstrate that the EnKF-produced w is much better than that produced by 4D-Var.To understand why a good EnKF-producedw leads a bad c
1
w, we examine an experiment for details with the observational interval, ensemble size, and MTW being 5 time steps, 50 time steps, and 500 time steps. Figure 11 plots the time series of w (Figure 11(a)) and c
1
w (Figure 11(b)) for 4D-Var CPE (black), truth (red), and EnKF CPE (blue) in EXP-1. Obviously, the EnKF-produced w is better than 4D-Var-produced w. We start from the following formulation of the absolute error of c
1
w to answer the question at the beginning of this paragraph:(12)
c
1
e
s
t
w
e
s
t
-
c
1
t
r
u
w
t
r
u
,where the superscripts “est” and “tru” indicate the estimation and truth of c
1 and w. As we stated in Section 2.2, the truth value of c
1 is 0.1 while the initial biased value of c
1 in 4D-Var CPE is 0.11 and ensemble means of c
1 for all ensemble sizes in EnKF CPE also approximate 0.11. Therefore, the computational process of (12) can be described as first scaling the truth value of w by 0.1 and the values of estimated w values by 0.11, then computing the absolute value of the difference. From Figure 11, it happens that the first step converts the worse 4D-Var-produced w to the appearing better c
1
w. Therefore, it is difficult to determine which algorithm produces better coupling processes without parameter estimation.Figure 11
Time series ofw (panel (a)) and c
1
w (panel (b)) for 4D-Var CPE (black) with MTW being 500 time steps, truth (red), and EnKF CPE (blue) with ensemble size being 50 in EXP-4. Here, the observational interval is set to 5 time steps.
(a)
(b)To fairly compare the performances of two algorithms, we focus on the results in EXP-4. From Figures4(d) and 8(d), we find that the representation of the sea-to-air coupling process produced by 4D-Var CPE with a MTW over 500 time steps is comparable to EnKF CPE. However, the EnKF CPE with ensemble size less than 50 improves the assimilation quality of 4D-Var CPE with a MTW shorter than 500 time steps by 59% (from 0.157 to 0.065) on average. Here, we choose an extreme case (i.e., the observational interval, ensemble size, and MTW being 50 time steps, 10 time steps, and 50 time steps) to compare the performance of two algorithms. Figure 12 shows time series of absolute errors of c
1
w for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. The EnKF CPE is better than the 4D-Var CPE. According to Figures 3(d) and 7(d), in the above parameter setting, the EnKF-estimated c
1 is much better than the 4D-Var-estimated c
1, since the MTW is too short to exactly retrieve c
1 for 4D-Var. In addition, we found that the 4D-Var-produced w is comparable to that produced by EnKF (not shown). Thus, the superiority of EnKF CPE over 4D-Var CPE is mainly attributed to the estimation of c
1.Figure 12
Time series of absolute errors ofc
1
w for 4D-Var CPE (blue) with MTW being 50 time steps and EnKF CPE (red) with 10-ensemble size in EXP-4. Here, the observational interval is 50 time steps.For the air-to-sea coupling process, comparison between Figure5(d) and Figure 9(d) reveals that 4D-Var CPE outperforms EnKF CPE when MTW exceeds a critical value (here is about 500 time steps) while EnKF CPE has advantages over 4D-Var CPE with a short MTW for a small observational interval. Here, we take two extreme cases to compare the performances of 4D-Var CPE and EnKF CPE. First, we check the results of two methods for the observational interval, ensemble size, and MTW being 50 time steps, 100 time steps, and 1000 time steps. Figure 13(a) plots time series of absolute errors of c
2
x
2 for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. After the spin-up period (about 40 TUs), the air-to-sea coupling process represented by 4D-Var CPE is better than that produced by EnKF CPE. Although the observational interval is very large and both 4D-Var CPE and EnKF CPE can effectively estimate c
2 (see Figures 3(b) and 7(b)), 4D-Var CPE can produce better x
2 than EnKF CPE (not shown). Second, we check the results (Figure 13(b)) of two methods for the observational interval, ensemble size, and MTW being 5 time steps, 10 time steps, and 50 time steps. According to Figures 3(b) and 7(b), it is difficult for 4D-Var CPE to estimate c
2 with a short MTW while EnKF CPE can well estimate c
2 with an ensemble size of 10. For x
2, 4D-Var CPE can produce similar results as EnKF CPE (not shown). Therefore, the EnKF CPE-produced c
2
x
2 is better than that produced by 4D-Var CPE.Figure 13
Time series of absolute errors ofc
2
x
2 for 4D-Var CPE (blue) and EnKF CPE (red) in EXP-4. In panel (a) (panel (b)), values of the observational interval, ensemble size, and MTW are 50 (5) time steps, 100 (10) time steps, and 1000 (50) time steps.
(a)
(b)


## 6. Summary and Discussions

A simple coupled model that characterizes the interaction of media with two different time scales is used to study the feasibility of the 4D-Var and EnKF coupling parameter estimation (CPE). Within a perfect OSSE framework which assumes that model errors only arise from the erroneously-set coupling parameters, the results demonstrate that, compared to traditional state estimation, both 4D-Var CPE and EnKF CPE algorithms can greatly improve the representation of air-sea coupling processes. An appropriate MTW exists in the 4D-Var CPE. Thus, if 4D-Var is used to implement CPE, the cost function should include observations from both the atmosphere and ocean, and an appropriate MTW should also be chosen. A roughly comparison between two assimilation methods demonstrates that EnKF CPE outperforms 4D-Var CPE with a MTW shorter than 5 TUs (i.e., 180 model days) for the sea-to-air coupling process. For the air-to-sea coupling process, 4D-Var CPE with a MTW longer than 5 TUs is better than EnKF CPE. The EnKF CPE is better than the 4D-Var CPE with a MTW shorter than 5 TUs for short observational intervals.Many challenges remain before CPE can be applied to a coupled general circulation model (CGCM). First, the perfect OSSE that uses the same model to setup the Truth (also called the Nature Run) and the assimilation experiments can be overoptimistic in the results. Thus, the perfect OSSE should be first extended to a biased OSSE and then put forward to the real world. Second, some of the above conclusions depend on the sensitivities of the coupling parameters and model errors. Thus, the conclusions should be further validated under more complicated model errors in CGCMs. Third, the static-multiplicative inflation scheme in the EnKF CPE should be updated to an adaptive inflation scheme [44]. Last, the impact of the observing system on the CPE should be examined with a GCM or a CGCM through varying the temporal and spatial densities of different observation variables or types.



---

*Source: 530764-2015-08-06.xml*