## Abstract

Imperfect physical parameterization schemes in a coupled climate model are an important source of model biases that adversely impact climate prediction. However, how observational information should be used to optimize physical parameterizations through parameter estimation has not been fully studied. Using an intermediate coupled ocean–atmosphere model, the authors investigate parameter optimization when the assimilation model contains biased physics within a biased assimilation experiment framework. Here, the biased physics is induced by using different outgoing longwave radiation schemes in the assimilation model and the “truth” model that is used to generate simulated observations. While the stochastic physics, implemented by initially perturbing the physical parameters, can significantly enhance the ensemble spread and improve the representation of the model ensemble, the parameter estimation is able to mitigate the model biases induced by the biased physics. Furthermore, better results for climate estimation and prediction can be obtained when only the most influential physical parameters are optimized and allowed to vary geographically. In addition, the parameter optimization with the biased model physics improves the performance of the climate estimation and prediction in the deep ocean significantly, even if there is no direct observational constraint on the low-frequency component of the state variables. These results provide some insight into decadal predictions in a coupled ocean–atmosphere general circulation model that includes imperfect physical schemes that are initialized from the climate observing system.

## 1. Introduction

Coupled climate model results are usually biased from those from the real world due to the imperfect numerical implementations and physical parameterization schemes as well as improper parameter values (Zhang et al. 2012). These model biases maintain systematic errors in the traditional climate estimation and produce artifacts in the analyzed variability (e.g. Dee and Da Silva 1998; Dee 2005). The model prediction in a biased model with traditional initialization schemes tends to drift away from observed states and thus has limited forecast skill (Smith et al. 2007).

Regardless of whether the model error arises from dynamical core misfitting, physical scheme approximation, or model parameter errors, the model errors, in principle, could be corrected using data assimilation. Using observations to optimize uncertain parameters in climate models is an important and complex subject. Many efforts have been undertaken to include model parameters into data assimilation control variables (e.g. Daley 1991; Wunsch 1996; Anderson 2001; Bennett 2002; Annan et al. 2005; Aksoy et al. 2006; Evensen 2007; Kondrashov et al. 2008) using four-dimensional variation or ensemble Kalman filter data assimilation methods. Some studies have addressed parameter optimization in a coupled climate model in which the interactions of multiple time scales play important roles in the development and propagation of climate signals. For example, a coupled data assimilation scheme with enhancive parameter correction (DAEPC) is designed to determine how to obtain a signal-dominant state-parameter covariance in order to effectively optimize the coupled model parameters using observations in different system components (Zhang et al., 2012).

Coupled model parameter estimation that uses observational information to tune the coupled model parameters has shown great potential to reduce model biases and improve the quality of climate estimation (i.e., model state estimation) and prediction (i.e., model predictability; Zhang, 2011a,b). Using a simple pycnocline prediction model, Zhang (2011a) investigated the impact of coupled state-parameter optimization on decadal-scale predictions, and showed that when both the coupled model states and parameters are optimized with data assimilation, the coupled model’s predictability is greatly improved. With a DAEPC algorithm in an intermediate coupled model, Wu et al. (2012) introduced a geography-dependent parameter optimization (GPO) scheme to increase the signal-to-noise ratio in the parameter estimation in a perfect model twin experiment framework, and examined the impact of the scheme on climate estimation and prediction (Wu et al. 2013). Recently, the DAEPC algorithm has also been demonstrated successfully in a coupled ocean–atmosphere general circulation model (Liu et al. 2014).

However, all the coupled model studies mentioned above applied the DAEPC algorithm in a perfect model framework. As a result, the impacts of imperfect physical parameterization schemes, which are an important source of model biases, have not been examined. In a coupled model, physical processes are usually approximated by parameterization schemes, which could be biased in their structure due to incomplete understanding and/or representation of the physical processes. In this study, we try to answer this question: To what degree can climate estimation and prediction be improved through physical parameter optimization using observations to minimize the error induced by a biased model structure? As a pilot study, we first construct an intermediate coupled model, which includes physical parameterizations of the heat exchanges among atmosphere, land, and ocean. With the intermediate coupled model with biased physics, a biased assimilation experiment framework is designed to study parameter optimization using a DAEPC algorithm. The biased physics in this study is defined using different outgoing planetary longwave radiation schemes in the assimilation and truth models. Simulated observations sampled from the truth model are assimilated into the assimilation model, and the degree to which the assimilation and prediction, with or without physical parameter optimization, recovers the truth is a measure of the impact of physical parameter optimization on climate estimation and prediction.

Based on such an experimental format, we compare results from five different assimilation experiments: 1) state estimation only (SEO) in which only observations are assimilated into the model state; 2) SEO with perturbed parameters (SEO_PP) in the biased physical optimization scheme; 3) single-valued parameter estimation (SPE) in which each parameter in the biased physical scheme is optimized with a globally uniform value (PP_PO1); 4) SPE in which only the parameters in the biased physical scheme to which the state variables are most sensitive are optimized (PP_PO2); and 5) geography-dependent parameter optimization (i.e., GPO) in which parameters to which the state variables are most sensitive are optimized according to local observational information and model sensitivity thus allowing the parameters to vary geographically.

The paper is organized as follows. The formulation of the intermediate coupled model is detailed in section 2. Section 3 discusses the format of the biased model framework. Section 4 presents the development of model biases due to the use of different outgoing planetary longwave radiation schemes. The impacts of parameter optimization in the biased physical scheme on climate estimation and prediction are examined and discussed in sections 5 and 6, respectively. Finally, a summary and discussion are provided in section 7.

## 2. An intermediate coupled model

As a first step in addressing the parameter optimization in an assimilation–prediction model that includes a biased physical scheme, we construct an energy-conserving intermediate coupled ocean–atmosphere model. The model is a combination of the intermediate coupled model without physical parameterizations presented by Wu et al. (2012) and the globally resolved energy balance (GREB) model used by Dommenget and Flöter (2011). It should be noted that despite some limitations in the representations of basic processes in the GREB model, such as the absence of an ENSO dynamical mechanism, the GREB model can reasonably simulate aspects of the Arctic winter amplification, the equilibrium land–sea warming contrast, and the interhemispheric warming gradient in the equilibrium 2×CO_{2} response. In addition, the surface layer humidity equation in the GREB model is removed to not unduly complicate the coupled model constructed for this study.

The intermediate coupled model we construct includes a vorticity advection equation that represents the atmosphere, a 1.5-layer baroclinic ocean with two simple temperature equations for the surface (Liu 1993) and subsurface ocean, and a simple land surface temperature model to complete the bottom boundary condition for the atmosphere over the entire globe. The GREB model is based on a surface energy balance through simple representations of solar and thermal radiation as well as other associated physical processes. The detailed description and the source code of the GREB model can be found online at http://users.monash.edu.au/~dietmard/content/GREB/.

The complete intermediate model includes six prognostic variables to describe the evolution of the three model domains: atmosphere, land, and ocean. The atmospheric variables include the geostrophic atmospheric streamfunction, , and the temperature of the troposphere above the boundary layer, . The land is only described by the surface temperature,, which represents the 2-m air temperature over land. The oceanic variables include the oceanic streamfunction , the sea surface temperature , and the deep ocean temperature . The model equations are briefly described in the following subsections.

### a. The atmosphere

The atmosphere consists of a global, barotropic, spectral model based on an equation of potential vorticity conservation and an atmospheric temperature equation:

In the global barotropic spectral model equation (1a), , , is the Coriolis parameter, *y* denotes the northward (meridional) distance from the equator, is the surface flux coefficient, and is the conversion coefficient from atmospheric streamfunction to temperature.

In the atmospheric temperature (*T*_{a}) equation (1b), is the atmospheric heat capacity, denotes the isotropic diffusion coefficient in the atmosphere, is the wind field (the eastern and northern components of which are defined as and , respectively), and is the sensible heat flux defined by

where is a coupling constant and is either (*T* over land) or (*T* over ocean) depending on whether the grid points are over land or ocean, respectively. Also, is the latent heat released by condensation of atmospheric water vapor and given by

where is the constant latent heat of evaporation of water, and denotes the change per unit time in due to precipitation and given by

where = −0.1/24 h and is the humidity of the surface layer derived from the climatologic NCEP mean from 1950 to 2008. Furthermore, is the net upward longwave radiative flux into the atmosphere, including the longwave radiation emitted by Earth’s surface that is absorbed by the atmosphere and reemitted downward and defined as

where is the effective emissivity (which will be discussed in section 3b), = 5.67 × 10^{−8} W m^{−2} K^{−4} is the Stefan–Boltzmann constant, and is the atmospheric radiative temperature as defined in section 3.4 of Dommenget and Flöter (2011). Finally, is the outgoing planetary longwave radiation, which will be discussed in section 3b.

### b. The ocean

The ocean is represented by equations for the upper oceanic streamfunction , the sea surface temperature *T*_{s}, and the deep ocean temperature *T*_{o}:

In the equation of the oceanic streamfunction (2a), is the oceanic deformation radius, denotes the coupling coefficient between the atmospheric and sea surface temperatures, and is the horizontal mixing coefficient of .

In the sea surface temperature equation (2b), is the oceanic heat capacity, is the oceanic current field (the eastern and northern components of which are given by and , respectively), is the horizontal diffusivity coefficient of , and represents the strength of upwelling (downwelling). The heat flux terms , , , and are the absorbed incoming solar radiation, the net thermal radiation, the latent heat release to the surface layer, and the sensible heat flux, respectively. These are defined as follows.

where and are the cloud and surface albedo, respectively, is the solar constant, and denotes the 24-h mean fraction crossing a normal surface area at the top of the atmosphere, which is assumed to be a function of latitude, , and time, *t*.

where the variables are those defined in section 2a. Note that was also been described in the previous section.

where is the air density, is the transfer coefficient, is the wind speed derived from the geostrophic atmosphere streamfunction , and is the saturation specific humidity of the surface air layer [see Eq. (6) in Dommenget and Flöter 2011]. Finally, is the net heat flux between the deeper ocean and the surface layer, which is defined as

where is the heat exchange between the surface and the subsurface ocean given by , where is the transfer coefficient and indicates the entrainment of subsurface water into the mixed layer represented by

where and are the mixed layer depth and its change, respectively; is the empirical flux correction term of .

The equation for the deep ocean (2c) describes the evolution of the subsurface temperature field, simulated by a simple linear equation forced by the upper ocean, where is the time step size and is the heat capacity of deep ocean. The term denotes the local entrainment of surface layer water, which is given by

where is assumed to be and indicates that the heat exchange with the deeper abyssal ocean is not simulated here.

### c. The land

The land in the coupled model is represented simply by an equation for the temperature 2 m above the land surface that is driven by heat exchanges with the atmosphere:

where is the heat capacity of the land surface; , , , and are as described in section 2b, except for the heat exchange between the atmosphere and the land, and is the empirical flux correction term for .

### d. Coupling of the model components

The coupling between model components is accomplished mainly through the physical parameterizations of the flux exchanges rather than only by coupling coefficients as in Wu et al. (2012). Dommenget and Flöter (2011) described the parameterizations in detail for solar radiation, thermal radiation, sensible heat flux, latent heat flux, and other terms associated with energy balance. Although the physical processes are greatly simplified relative to a state-of-the-art CGCM, the coupled model is sufficient for the purpose of studying the ensemble coupled parameter optimizations with biased physics.

For simulating a mean climate as close to the true climatology as possible, flux correction terms are imposed on the equations of and . These terms are estimated by calculating the residual tendencies needed to maintain the observed climatology. A potential disadvantage of the inclusion of the flux correction terms is that they may affect the parameter optimization in a nonlinear way. However, Dommenget and Flöter (2011) pointed out that many studies have addressed climate sensitivities using flux corrections in a simplified ocean model and most have found that the climate sensitivity is close to that without flux correction. We will discuss this issue further in section 5b. But it should be pointed out here that, despite our attempt to relate this model to the real world, the model here is meant to serve mainly as a conceptual model of a coupled system with different time scales and high dimensionality. Therefore, the results obtained from this model should have some general implications for more realistic coupled ocean–atmosphere models.

### e. Model mean states and variability

The basic configuration of the coupled model, such as the horizontal resolution, numerical time step, and the rate of change of the greenhouse gases, are consistent with those of the GREB model. The default values of all the parameters are set as in Table 1 of Wu et al. (2012) and Table 2 of Dommenget and Flöter (2011). Starting from the initial conditions,

the coupled model is integrated for 350 years, where *ψ*^{0}, , , , , and are the corresponding climatological fields of the seven prognostic variables. The annual means of *ψ*, , , , , and over the last 300 years of the simulation are shown in Figs. 1a–g (note that the first 50 years were considered as spinup time and discarded). The first thing to note is that these results contain a reasonable atmospheric temperature gradient and associated wave trains in the *ψ* field are observed. In the field, the large-scale ocean circulation is characterized by gyre circulations in the Pacific, Atlantic, and Indian Oceans as well as an Antarctica Circumpolar Current (ACC) in the Southern Ocean. Because of the coupling, the patterns of and are very similar to the atmospheric temperature over the land and ocean, respectively. Since is forced directly via the entrainment and vertical turbulent mixing between the sea surface and the subsurface, the pattern of is also similar to .

The annual-mean absolute errors of the sea surface temperature and the land surface temperature are calculated as the difference from the model results to their respective climatologic datasets. As can be seen from Fig. 2a, the absolute values of the error in in the most regions of the Pacific, Atlantic, and Indian Oceans are less than 5 K (see the white and red regions in Fig. 2a). The errors in the surface temperature show high values around the globe in the Arctic Ocean and the ACC region, with values in some locations exceeding 10 K. For , the absolute values of the error over Asia, North America, South America, and West Antarctica (0°–180°E) are less than 10 K (see the white and red domains in Fig 2b). The errors in remain large (exceeding 10 K) in some regions such as Africa, Australia, and eastern Antarctica (180°–360°E). Although some ocean and land regions present visible biases relative to the climatology, Dommenget and Flöter (2011) indicated that in many studies with simplified models modest mean state errors (less than 10 K) are acceptable and assumed to be of minor importance.

To investigate the variability of the state variables, an empirical orthogonal function (EOF) decomposition of the 300-yr time series of the monthly mean anomalies of *T*_{a} was conducted. The time coefficients of the first four EOFs were used to perform the power spectrum analysis to show the internal variability of *T*_{a}. The characteristic time scales of the first four modes were 100 years, 50 years, 10–25 years, and 6–14 years (Fig. 3), indicating that the model states consist of multiple time scales. Figure 4 shows the spatial distributions for the first four EOF modes of the monthly averaged atmospheric temperature . Note that Fig. 4 does not show an obvious ENSO cycle because this simple coupled model lacks the ability to simulate the variability of the thermocline depth in the ocean. However, despite limitations in the representations of some basic physical processes, the model is of sufficient mathematical complexity for the purposes of this study of parameter optimization with an assimilation model that contains biased physics.

## 3. Experimental setup

In this section, using the intermediate model described above and an ensemble coupled data assimilation scheme, we designed a biased model framework by setting outgoing planetary longwave radiation schemes different between the assimilation model and the truth model. We assume the different parameterizations of outgoing longwave radiation to be the only source of assimilation model biases aside from errors in the initial model states.

### a. Ensemble coupled data assimilation with enhancive parameter correction

A coupled data assimilation scheme with enhancive parameter correction (i.e., DAEPC) (Zhang et al. 2012) is employed to compute the model states and perform parameter optimization. This is a modification to the standard data assimilation with adaptive parameter estimation scheme described by Kulhavy (1993) and Tao (2003). The covariance between a parameter and the model state serves as the critical quantity to project observational information of the model states to the parameter being optimized.

Based on a two-step ensemble adjustment Kalman filter (EAKF; Anderson 2003; Zhang and Anderson 2003), the observational increment is computed first. Equation (4) indicates that the new ensemble spread, , is calculated by the prior ensemble spread, , and , where *i* represents the ensemble sample index and *k* represents the observation index; is determined by Eq. (5), where and represent the standard deviation of the observation and its prior estimate from the ensemble, respectively. The shift of ensemble mean induced at the *k*th observational position, , is computed according to Eq. (6), where and are the prior ensemble mean of the estimated value and observation value at the *k*th position, respectively. Then the increment induced by the observation, , is obtained using Eq. (7), in which is the prior estimated value at the *k*th observational position for the *i*th ensemble member.

After the observational increment, , is calculated, it is projected onto the corresponding model variables and model parameters using uniform linear regression formulas:

Note that Δ*x*_{j,i} indicates the contribution of the *k*th observation to the state variable, *x*, at the *j*th grid point for the *i*th ensemble member; indicates the contribution of the *k*th observation to the model parameters, , at the *j*th grid point for the *i*th ensemble member. If each parameter being estimated is assumed to be a globally uniform value at all grid points, and will degenerate to and , respectively. The term denotes the error covariance between the prior ensemble of the state variables (the parameters) and the model-estimated ensemble of the observation in the *k*th observational position. The inflation scheme used in DAEPC by Wu et al. (2012) is also adopted in this study. Additionally, the distance factor (Zhang et al. 2007) is introduced into the DAEPC scheme to remove spurious correlations over long distances. The impact radius of observations is set to 500 km for *ψ* and and is set to 1000 km × for and , where is the latitude of model grid points.

Unlike the model states, the model parameters in the coupled climate model do not have any dynamically generated internal variability. Therefore, the quality of the covariance between a parameter and the model state [see Eq. (9)] is the only criterion needed to determine if the parameter can be optimized reasonably well. The subsequent parameter updating with observational data can improve the state estimation of the next cycle, and this improved state estimation further enhances the signal-to-noise ratio of the state parameter covariance for the subsequent cycle of parameter correction.

### b. Two schemes of outgoing planetary longwave radiation

In our ocean–atmosphere–land coupled climate model, the outgoing planetary longwave radiation is the only heat loss from the coupled system. The magnitude of the heat loss determines how much heat remains in the entire system, and this process further modulates the distribution of the heat reserves within the coupled system. Here, we introduce two parameterization schemes for outgoing planetary longwave radiation. The first scheme considers the outgoing planetary longwave radiation, , of the planet if it acts as a gray body with emissivity . Fanning and Weaver (1996) also used this form of to develop an atmospheric energy–moisture balance model. Note that in the GREB model is given by

where is the effective emissivity, is the Stefan–Boltzmann constant equal to 5.67 × 10^{−8} W m^{−2} K^{−4}, and is the atmospheric radiative temperature defined in section 3.4 of Dommenget and Flöter (2011). The effective emissivity in Eq. (6) is given by

where CLD is the cloud cover, taken from NCEP reanalysis; is the emissivity without clouds given by

where CO_{2} is the atmospheric concentration of carbon dioxide and is the vertically integrated atmospheric water vapor calculated by Eq. (8) of Dommenget and Flöter (2011). The 10 parameters are taken from Table 2 of Dommenget and Flöter (2011), hereafter denoted the DF11 scheme.

The second scheme for outgoing planetary longwave radiation is taken from Thompson and Warren (1982) but modified to include the changes in CO_{2} concentration and formulated as

where *r* is the surface relative humidity calculated according to the addendum of Thompson and Warren (1982) and is the CO_{2} concentration at the reference level ( = 350 ppm). Note that *F* = 5.77 W m^{−2} corresponds to a specified radiative forcing of 4 W m^{−2} for a doubling of atmospheric CO_{2} (Ramanathan et al. 1987). The values of the parameters to can be found in Table 3 of Thompson and Warren (1982) and for convenience of discussion they are listed in Table 1 herein. Hereafter, we refer to this scheme as the TW82 scheme. The TW82 scheme was used by Weaver et al. (2001) to construct an Earth system climate model.

The two physical schemes, DF11 and TW82, have completely different implementations of the same physical processes. The TW82 scheme uses a cubic polynomial in to curve fit a detailed radiative transfer model, and the coefficients of the polynomial fit take quadratic forms with respect to surface relative humidity. By contrast, the DF11 scheme has a simple dependence on the atmospheric radiative temperature that treats the atmosphere as a grey body, but the emissivity coefficient is rather complicated, being determined by 10 empirical parameters and associated atmospheric variables.

### c. Data assimilation experiment with biased physical schemes

The truth model uses the DF11 scheme. Starting from the initial conditions described in section 2e, the truth model is run for 100 years to generate the time series of the model truth with the first 50-yr considered the spinup period. The model observations of , , , and are generated by sampling the state of the truth model at a specific observational interval. Gaussian white noise is added to the model observation to simulate observational error. The standard deviations of the Gaussian observational errors are 0.5 K for and , 1 K for , and 10^{6} m^{2} s^{−1} for . The sampling intervals are 12 h for , , and and 48 h for . All observation locations are randomly distributed globally with the same density as the model grid.

The assimilation model uses the TW82 scheme. Starting from the same initial conditions used in the truth model run, the assimilation model is spun up for 50 years to produce the biased initial conditions . Based upon , the ensemble initial conditions are created by adding Gaussian white noise to with a standard deviation of 10^{6} m^{2} s^{−1}. Consistent with previous studies (Zhang 2011a,b; Wu et al. 2012, 2013), the ensemble size is set to 20 and we denote the ensemble initial conditions as .

Starting from the ensemble initial conditions, , five assimilation experiments are conducted. The details of all five assimilation experiments and model simulations are listed in Table 2. First, the experiment for the state estimation only, which only assimilates observations into model states, is performed. Second, the SEO with the 12 parameters listed in the Table 1 initially perturbed with a standard deviation of 5% of the default values is performed and denoted SEO_PP. Third, the single-valued parameter estimation is performed to optimize each parameter listed in the Table 1 with globally uniform values, which is denoted PP_PO1. Fourth, the SPE is performed to only optimize parameters to which is the most sensitive, denoted PP_PO2; Finally, the geography-dependent parameter optimization is carried out to optimize the parameters to which is the most sensitive. The inflation scheme used in GPO follows that used in Wu et al. (2012). In addition, a control run without any observational constraints, denoted CTRL, is performed to serve as a reference for the evaluation of the assimilation experiments.

It should be pointed out that the assimilation model differs from the truth model only in the longwave scheme employed [Eq. (13) vs. Eqs. (10)–(12)]. Because of the different parameterization schemes for the longwave radiation, there is no longer a true value for the parameters in the longwave scheme in the assimilation model Eq. (13). The only metrics to evaluate the performance of a parameter estimation scheme are the accuracy of the estimated model states and/or the forecast. This configuration simulates a real world scenario, where the parameters in a climate model have no true values, so that parameter optimization can only be evaluated by examining the accuracy of the estimation and model prediction.

## 4. The model bias arising from physical schemes

This section examines the model biases induced by different physical schemes, that is, the difference between the annual mean and climatology in the assimilation and truth models. Figs. 5a–d show the spatial distribution of biases in the annual mean of *ψ* , , , and . First, is directly affected by the different parameterization of the outgoing planetary longwave radiation in the two schemes. An obvious mean bias in can be seen in Fig. 5b, which manifests as a cold bias in high latitudes and over land in the northern middle latitudes (including the Eurasian continent and the Northern American continent). The maximum value of this cold bias exceeds 10 K and is seen to be larger in the Northern Hemisphere (NH) than the Southern Hemisphere (SH). In the SH, the largest biases occur over land, whereas in the NH they occur over the Arctic Ocean. This latter bias may be associated with the absence of sea ice dynamical process. The biases in over land are larger than over the ocean, indicating that over land is more sensitive to the parameterization of outgoing planetary longwave radiation than that over the ocean. The bias in in the tropical ocean and land is small (<2 K), and the largest warm bias lies in the southeast region of Australia and is more than 5 K.

Note that affects through the interaction between the atmosphere and ocean, and therefore the pattern of the bias in is similar to that in (see Fig. 5c). Further, since the deep ocean temperature is controlled by , via the local entrainment and the vertical turbulent mixing between the sea surface and the deeper ocean, the distribution of the bias in (see Fig. 5d) is similar to that in .

The bias in *ψ* (see Fig. 5a) is strongly negative around the globe in the polar and subpolar regions, especially in the NH high latitudes. In the low and middle latitudes, the bias in *ψ*, while still mostly negative, is weaker and the patterns are more complex. This complexity is partly due to the stronger internal variability of *ψ* in the tropical and subtropical regions.

Figure 6 shows the climatologic annual cycle of global mean *ψ*, , , and in the truth model run (dashed line) and the assimilation model run (solid line). The mean biases in *ψ*, , , and are about 2 × 10^{5} m^{2} s^{−1}, −7 K, −3 K, and −2 K, respectively. The of the assimilation model is lower than that of the truth model, indicating that the TW82 parameterization scheme in the assimilation model [Eq. (13)] leads to more longwave radiation being emitted to space than that from the DF11 scheme in the truth model [Eqs. (10)–(12)]. This also leads a lower spatial averaged and in the assimilation model.

## 5. Parameter optimization with biased model physics

The analyses in the previous section showed that model climate depends heavily on the representation of model physics, even when the flux correction terms work correctly in the coupled model. Differences in outgoing planetary longwave radiation schemes create different model climatology. In this section, we will examine how this difference affects climate estimation and prediction, and how to overcome this problem by optimizing the parameters in TW82 scheme using observations produced by the DF11 scheme. We begin by examining the sensitivities in the assimilation model with respect to the parameters in the TW82 scheme. Then, based upon the model sensitivity, we discuss the impact of parameter optimization on climate estimation.

### a. Model sensitivities with respect to physical parameters

It is essential to investigate model sensitivities with respect to parameters prior to parameter optimization. In this study, the sensitivity study is carried out for all 12 parameters in the TW82 scheme (Table 1). The ensemble spread of a model prognostic variable is used to quantitatively evaluate the relative sensitivities. For each parameter, 20 random numbers are drawn from a normally distributed population with a standard deviation of 5% of the default value. This random number is superimposed on the default value of the parameter being perturbed, while other 11 parameters remain fixed at their default values. All 20 ensemble model runs are started from the same initial conditions, *S*_{1}, and the assimilation model is integrated for 70 years. Sensitivities are calculated using the model output from the last 50 years. This process is repeated for each parameter.

Figures 7a–f show the time-space averaged sensitivities of *ψ*, , , , , and with respect to all 12 parameters. For (Fig. 7c), the sensitivity with respect to all 12 parameters is more than 1.5 K, indicating that a slight perturbation in any parameter can produce a significantly different model state. We see that is more sensitive to the first parameter (*b*_{00}) and the last four parameters (*b*_{20}*–b*_{23}) than the other parameters. The sensitivity of to each of these five parameters is more than 2 K. Because of the strong interaction between the atmosphere and both the land and the ocean, parameter sensitivities of and are similar to those of . Further, due to the strong impact of on , the sensitivities of are similar to those of . Because of strong nonlinearity and internal variability of flow fields, the sensitivities of *ψ* and show some differences from the sensitivity of .

The parameter sensitivities of the model states vary significantly with geographic location. For example, Fig. 8a presents the distribution of time-averaged sensitivity of to the parameter *b*_{00.} The most sensitive areas of are over the land regions in the NH high latitudes, eastern Asia, the eastern part of the North Pacific Ocean, and in the Antarctic Circumpolar Current region. This indicates that a slight perturbation in *b*_{00} can induce a dramatic drift of in those regions. In contrast, is largely insensitive to perturbations in *b*_{00} in the tropical regions, the southern Africa and over the continent of Australia. The geographic patterns of the sensitivity of the model state to the other 11 parameters (not shown) have patterns similar to that shown in Fig. 8a.

The relative sensitivity, defined as the sensitivity normalized by the standard deviation, is also an interesting quantity that explains the significance of the sensitivity of a model variable. For example, Fig. 8b shows the distribution of time-averaged relative sensitivity of . Contrary to the sensitivities shown in Fig. 8a, the relative sensitivity over the land regions in high latitudes is lower than 0.2, which indicates that has a strong internal variability over this region. In contrast, the relative sensitivity of over the ocean is larger than that over land, especially in the tropics and south of Australia, where the values all exceed 1.0, suggesting that the internal variability of over ocean is smaller than that over land.

### b. Impact of parameter optimization on climate estimation

Figures 9a–d show time series of the RMSEs of *ψ*, , , and for the CTRL (pink line), SEO (blue line), SEO_PP (red line), and PP_PO1 (black line) experiments. It should be noted that all RMSEs calculated in this study are priors. From Fig. 9b, we can see that the average RMSE of for the CTRL experiment reaches about 9 K, with an annual cycle amplitude of ~2 K, which shows the obvious difference induced by the two different longwave radiation schemes. The variability of in the SEO experiment shows significant improvement (blue line in Fig. 9b) when compared to CTRL, in that the average RMSE decreases to ~4 K, with an annual cycle amplitude of about 0.5 K. The average RMSE of in the SEO_PP experiment (red line in Fig. 9b) is reduced further by 59% compared to that of SEO, indicating that the stochastic physics through perturbing the physical parameters initially can significantly improve the quality of the state estimation. Recall from section 5a that the perturbation of the parameters can increase the spread of model states, which increases the representation of the ensemble for the model uncertainty. Relative to SEO_PP, the PP_PO1 experiment (black line in Fig. 9b) shows further reduction in the RMSE of . The RMSE of in PP_PO1 is reduced by about 55% relative to that of SEO_PP. We see that observational information is effectively extracted by the data assimilation process, and the model biases induced by the different physical parameterizations are efficiently reduced by observation-optimized model parameters.

The RMSE of in SEO is also significantly improved relative to CTRL (pink vs blue lines in Fig. 9c). The RMSE decreases from roughly 7 K to about 3 K, with the amplitude of the annual cycle cut roughly in half from 2 K to about 1 K. The amplitude of the oscillation in is larger than that in because the sampling frequency for is lower than that for , so that fewer observations of are assimilated into the model. The further improvement of in SEO_PP (the red line in Fig. 9c) shows that SEO_PP performs better than SEO. From SEO_PP, the parameter optimization experiment (PP_PO1) further reduces the RMSE by about 48%. Note that in PP_PO1, all 12 parameters are optimized from the observations of , so that the optimized parameters are modulating indirectly via the coupling of the atmosphere and the ocean.

Unlike the RMSEs of sea surface temperature and air temperature, which are reduced significantly during the first 100 days of assimilation, the RMSEs of in SEO, SEO_PP, and PP_PO1 decrease much more slowly, finally reaching equilibrium after 6000 days (Fig. 9d). This slow decrease of errors occurs because there is no direct observational constraint on . The RMSE of for PP_PO1 is the smallest of the three assimilation experiments, indicating the strong impact of parameter optimization on the deep ocean.

The RMSEs of *ψ* for the SEO, SEO_PP, and PP_PO1 experiments are all reduced by roughly two orders of magnitude relative to that for CTRL (see the small inset panel in Fig. 9a). However, the RMSE of *ψ* behaves differently than that for temperature in that the RMSE for PP_PO1 is not always smaller than that for the SEO or SEO_PP (blue, red, and black lines in Fig. 8a). There are two reasons for this. First, the RMSE of *ψ* for SEO has dropped to about 7 × 10^{5} m^{2} s^{−1}, which is already less than the standard deviation of observational errors of *ψ* (10^{6} m^{2} s^{−1}). Therefore, it is difficult to extract an observational signal to further reduce the RMSE of *ψ* through the adjustment of parameters*.* Second, physically, all 12 perturbed parameters in the longwave scheme directly impact the temperature field but not the streamfunction, *ψ*.

Spatial distributions of the RMSEs of in the SEO, SEO_PP, and PP_PO1 experiments are presented in Figs. 10a–c. The RMSE of over land is much larger than over ocean, especially in northern Africa, Asia, and Australia. This suggests that the temperature over land is affected much more by the biased physics than the temperature over the ocean. The SEO_PP experiment greatly reduces the RMSE of over land and over the Southern Ocean relative to SEO. Furthermore, the RMSE of in PP_PO1 is greatly reduced from SEO_PP in the most regions, especially over land. In general, the spatial distribution of RMSEs of in PP_PO1 depends on the spatial distribution of the sensitivity of . For example, the sensitivity of in western Antarctica is small (see Fig. 8a), so the reduction of RMSE is not significant there. In contrast, in eastern Asia has a large sensitivity, so the optimization of parameters there contributes significantly to the improvement of .

To further improve the performance of the parameter optimization, we performed the two additional assimilation experiments, PP_PO2 and GPO. Figures 11a–d show the time series of the RMSEs of , , , and in the last 10 years of 50-yr assimilation runs from PP_PO1 (pink line), PP_PO2 (blue line), and GPO (red line). The RMSE of has a seasonal cycle, which can be clearly seen from Fig 11a. The reduction of the RMSE is distinct, especially in boreal summer, suggesting that the contribution of the parameter optimization to the improvement of the model state is larger at that time of the year. We can see from the blue and pink lines in Figs. 11a–d that the improvement of the state estimation is greater when only optimizing the five parameters to which the state variables are most sensitive (PP_PO2) compared to the case of all parameters being optimized (PP_PO1). This suggests that the noise of parameter-state covariance in PP_PO1 is larger than that in PP_PO2, and that optimizing the parameters to which the state variables are most sensitive can better incorporate observational information into the model dynamics. The RMSE of in PP_PO2 also shows an obvious improvement compared to PP_PO1, especially in the boreal winter (blue vs pink lines in Fig. 11b). The RMSEs of and for PP_PO2 are also decreased relative to the RMSE in PP_PO1 (see Figs. 11c,d), showing a better quality of the parameter-state covariances over the ocean in PP_PO2.

The benefit of GPO is that it allows the optimized parameter values to vary spatially. This may be justified in that the difference between the assimilation model and truth may depend on the state climatology itself, which varies spatially, so that a spatially uniform parameter may not be the optimal representation of the physics. Compared to PP_PO2, GPO further reduces the RMSEs of , , , and (red lines in Figs. 11a–d). For example, the RMSEs of and in boreal winter are reduced by more than 20% and 50% respectively in GPO relative to that in PP_PO2, showing the great advantage of geography-dependent parameter optimization. The reduction of errors of is less significant than the reduction of errors, suggesting that is more sensitive to the parameters being optimized. Table 3 shows the total RMSEs of , , , and for all experiments. From the first column, PP_PO1, PP_PO2, and GPO, with reductions of about 56%, 61%, and 68%, respectively, all greatly reduce the RMS errors of relative to SEO_PP. Recall that SEO_PP already had dramatically reduced RMS errors in from SEO (by 59%). Further, GPO is also the best methodology to improve the estimation of , in that the RMSE of in GPO is reduced by about 17%, 31%, and 70% over those in PP_PO2, PP_PO1, and SEO_PP, respectively.

A perfect assimilation experiment without any model bias, SEO_PRFT, is also carried out to serve as a benchmark. SEO_PRFT is the same as SEO except that the assimilation model uses the DF11 scheme and all parameters in the DF11 scheme use the truth values. As can be seen from the solid black lines in Figs. 11a–d, the RMSEs of all four variables, , , , and for SEO_PRFT are the smaller than all of the other assimilation experiments. Of the three schemes PP_PO1, PP_PO2, and GPO, the RMSEs of , , , and for GPO are the closest to those for SEO_PRFT (red vs black lines in Figs. 11a–d), although there is still a visible difference. The difference indicates that the model bias due to the substantial difference in the two longwave parameterization schemes can be mitigated to some degree by parameter optimization.

To investigate the impact of the flux adjustment in the coupled climate model on the parameter optimization within our biased assimilation experimental framework, another four assimilation experiments are conducted. These are the same as SEO_PP, PP_PO1, PP_PO2, and GPO, but without the flux correction terms. Figure 12 shows the time series of the RMSEs of (Fig. 12a), (Fig. 12b), (Fig. 12c), and (Fig. 12d) for SEO_PP (black line), PP_PO1 (pink line), PP_PO2 (blue line), and GPO (red line) without flux corrections. In general, these RMSEs time series have similar behavior to those with flux adjustment (shown in Figs. 10a–d). That is, the RMSE decreases from SEO_PP to PP_PO1 to PP_PO2 and finally to GPO. Furthermore, the improvement in GPO is especially significant for the low-frequency signals in . Overall, this suggests that the flux adjustment in the coupled climate model has no significant effect on the ensemble coupled parameter optimization. Nevertheless, further studies are needed to better understand the effect of the flux correction on parameter optimization.

The GPO-optimized parameters have distinctive spatial patterns. For example, Fig. 13 shows the spatial distribution of *b*_{00}. The pattern shows that *b*_{00} is not far away from the default value of *b*_{00} (243.414) near the equator, but decreases towards high latitudes. In addition, *b*_{00} tends to be smaller over land than over the ocean (e.g. the Asia vs the Pacific Ocean and the Africa vs the Indian Ocean in the same latitude). Because of the functional form of the parameterization scheme, this means that *b*_{00} has a greater impact on the land than on the ocean, which is consistent with the greater reduction of RMSEs over land than over ocean discussed earlier. In addition, the biased physical scheme can cause larger model errors over the land and at high latitude than over the ocean or at low latitude. This is also consistent with the greater error reduction seen over land in the GPO scheme discussed earlier.

## 6. Impact of observation-optimized physical parameters on climate predictions

Although justified to some extent, the parameter optimization, as judged from the analysis, is subject to the criticism of curve fitting. However, from a more practical point of view, the role of parameter optimization should be judged from its effectiveness in climate prediction. In order to evaluate the impact of ensemble coupled data assimilation with biased model physics on model prediction, a 20-member ensemble forecast is performed. This forecast uses 20 initial conditions selected every 2 years apart from the SEO, SEO_PP, PP_PO1, PP_PO2, and GPO analysis fields over the last 40 years. The 20-member ensemble is integrated for 10 years starting from the analyzed states of each of the five assimilation schemes. The global mean of the anomaly correlation coefficient of the forecasted ensemble mean of and global RMSE is used to be the forecast skill evaluation metric. Further, we use to represent climate variability and to represent weather variability.

We examine the climate prediction skill question first. Figure 14 shows the variation of the ACC (Fig. 14a) and the RMSE (Fig. 14b) as a function of the forecast lead time for the forecasted ensemble means of for the SEO (green line), SEO_PP (black line), PP_PO1 (pink line), PP_PO2 (blue line), and GPO (red line) experiments. The initial ACC of for all five forecasted experiments are very high and close to 1.0. The ACC for SEO, SEO_PP, and PP_PO1 all drop below 0.75 after 120 months, but remain over 0.85 for both PP_PO2 and GPO. In addition, the RMSE of the GPO scheme is the smallest throughout the forecast lead time (Fig. 14b), suggesting that GPO has the best climate prediction skill. It is interesting that SEO_PP loses its superiority relative to SEO after about 80 months even though the initial state of SEO_PP is better than that of SEO. This may be caused by the noise in some of the parameters. Therefore, caution should be exercised when processing a long-term climate prediction using initially perturbed parameters in biased physical schemes without parameter optimization. Comparing SEO and GPO for the 120-month forecast lead time shows that the averaged ACC of in the GPO experiment is about 32% greater and the averaged RMSE is about 59% lower than those in SEO experiment.

As with the climate forecast, the parameter optimization also improves the weather forecast skill for the atmosphere. Figs. 14c and 14d show the variations of ACC (Fig. 14c) and RMSE (Fig. 14d) as a function of the forecast lead time for the forecast ensemble means of for SEO (green line), SEO_PP (black line), PP_PO1 (pink line), PP_PO2 (blue line), and GPO (red line). These figures indicate that the SEO scheme has the worst forecast skill, with the ACC dropping below 0.2 after 10 days into the forecast. SEO_PP increases the ACC and reduces the forecast error of the SEO scheme to some degree. GPO and PP_PO2 produce the best forecast performance after the 10-day forecast lead time, as they result in the highest ACC (~0.5) and the smallest RMSE (~4.5). Overall, the relative improvement in forecast skill delivered by parameter optimization is similar in both the weather and climate forecasts.

## 7. Summary and discussion

An imperfect dynamical core and empirical physical schemes and improper parameter values are three main sources of model bias. Among these sources of model error, in principle, only the model parameter errors can be treated directly by assimilating observations. Because of our incomplete understanding of atmosphere and ocean physics, we still have a long way to go to adequately improve the estimation and prediction capability of a coupled model through directly improving physical parameterization schemes. In this paper, we studied several methods for using observational information to optimize parameters in an attempt to minimize errors caused by a biased model structure. We designed a biased assimilation experiment framework to study parameter optimization in an energy-conserving intermediate coupled model. Within this framework, the assimilation-prediction model is subject to biased physics from the truth model. The biased physics was implemented by using different outgoing longwave radiation schemes in the assimilation–prediction model and the truth model. A series of assimilation and prediction experiments were performed to investigate the degree to which observational information can be used to optimize the physical parameterization and its associated parameters so as to improve climate estimation and prediction. While stochastic physics, implemented by initially perturbing the parameters of the physical parameterization scheme, can significantly enhance the ensemble spread and improve the representation of the model ensemble, the parameter estimation is better able to reduce model biases induced by the introduction of biased physics. Further, when the biased physics dominates the behavior of the coupled model, the parameter optimization where only the most influential parameters are optimized and allowed to vary geographically leads to better results for climate estimation and prediction. It is also noted that even when there is no direct observational constraint on the state variables of the low-frequency components, such as the deep sea temperature, the quality of the estimation/forecast of the state variables is improved with the DAEPC algorithm (especially the geography-dependent parameter optimization) compared to the traditional state estimation alone (SEO). In addition, compared to the state estimation within a perfect data assimilation framework, the DAEPC algorithm mitigates the model bias to some degree. Because of the substantial difference in the two physical parameterization schemes within the biased data assimilation framework, it is expected to improve the physical parameterizations themselves to mitigate the model bias further together with the DAEPC algorithm.

In spite of the promising results of parameter optimization in the intermediate climate model presented here, much more work is needed to fully understand the impact of model biases in coupled general circulation models (CGCMs) on real world climate estimation and prediction. Imperfect physical schemes are often used in CGCMs. Therefore, methods for optimizing model parameters in multiple biased physical schemes with the DAEPC algorithm need to be further examined to account for the possible compensating or noncompensating effects of the deficiencies in the physical schemes (Yang and Delsole 2009; Zhang 2011b), especially for the real world climate estimation and prediction.

## Acknowledgments

The authors would like to thank Drs. X. Yang, G. Vecchi, I. Held, Y.-S. Chang, and A. Wittenberg for their generous discussions. Suggestions from three anonymous reviewers contributed significantly to the final version of this work. This research is sponsored by NSF, 2012CB955200 and partly supported by grants from the National Basic Research Program of China (2013CB430304) and the National Natural Science Foundation of China (under Grants 41030854, 41106005, and 41206178).

## REFERENCES

*Inverse Modelling of the Ocean and Atmosphere.*Cambridge University Press, 234 pp.

*Atmospheric Data Analysis.*Cambridge University Press, 457 pp.

*Data Assimilation: The Ensemble Kalman Filter.*Springer Press, 279 pp.

*Adaptive Control Design and Analysis.*John Wiley & Sons, 640 pp.

*Climate Dyn.,*

**40,**1789–1798, doi:.

*The Ocean Circulation Inverse Problem.*Cambridge University Press, 442 pp.

**24,**6210–6226, doi:.