The retrospective prediction of El Niño-southern oscillation from 1881 to 2000 by a hybrid coupled model: (I) Sea surface temperature assimilation with ensemble Kalman filter

Abstract

In this study, we assimilated sea surface temperature (SST) data of the past 120 years into an oceanic general circulation model (OGCM) for El Niño-southern oscillation (ENSO) retrospective predictions using ensemble Kalman filter (EnKF). It was found that the ensemble covariance matrix in EnKF can act as a time-variant transfer operator to project the SST corrections onto the subsurface temperatures effectively when initial perturbations of ensemble were constructed using vertically coherent random fields. As such the increments of subsurface temperatures can be obtained via the transfer operator during assimilation cycles. The results show that the SST assimilation improves the model simulation skills significantly, not only for the SST anomalies over the whole assimilated domain, but also for the subsurface temperature anomalies of the upper 100 m over the tropical Pacific off the equator. Along the equator, the improvement of the assimilation is confined within the mixing layer because strong upwelling motions there prevent the downward transfer of SST information. The retrospective prediction skills of ENSO over the past 120 years from 1881 to 2000 were significantly improved by the SST assimilation at all leads of 1–12 months, especially for the 3–6 months leads, compared with those initialized by the control run without assimilation. The skilful predictions by the assimilation allow us to further study ENSO predictability using this coupled model.

Appendix: Estimation of the linear increment transfer operator

In this study, we consider the sea temperatures of the top 17 model levels. Thus X in (3) can be written X = [x ₁ x ₂ …x _k …x ₁₇ ]′, where x _k is the kth level’s temperature field. Denoted by m (m = 71 × 31)the number of spatial grids of assimilation domain at a level, x _k is an m-dimensional vector. P _e ^f in (3) is a covariance matrix between grids, which is an M × M(M = 17 × m) 2-dimentional matrix. We write P ^f _e using sub-matrices, the covariance between grids of two levels and the covariance between grids of the same level, namely,

$$ {{\mathbf{P}}}_{e}^{f}=\left[\begin{array}{ccccccc} {{\mathbf{C}}_{1,1} } & {{\mathbf{C}}_{1,2} } & \ldots & {{\mathbf{C}}_{1,k2}} & \ldots & {{\mathbf{C}}_{1,K-1} } & {{\mathbf{C}}_{1,K} }\\ {{\mathbf{C}}_{2,1} } & {{\mathbf{C}}_{2,2} } & \ldots & {{\mathbf{C}}_{2,k2}} & \ldots & \ldots & \ldots\\ \ldots & \ldots &&& \ldots\\ \ldots & \ldots &\ldots&\ldots& \ldots&&\ldots\\ {{\mathbf{C}}_{k1,1} } & {{\mathbf{C}}_{k1,2} } & \ldots & {{\mathbf{C}}_{k1,k2}} & \ldots & {{\mathbf{C}}_{k1,K-1} } & \ldots\\ \ldots & \ldots &&& \ldots&&\ldots\\ {{\mathbf{C}}_{K-1,1} } & \ldots& & \ldots & & {{\mathbf{C}}_{K-1,K}} & {{\mathbf{C}}_{K-1,K-1} } \\ {{\mathbf{C}}_{K,1} } & {{\mathbf{C}}_{K,1} } & & {{\mathbf{C}}_{K,j}} & & {{\mathbf{C}}_{K,K-1} } & {{\mathbf{C}}_{K,K} }\\ \end{array} \right]$$

(8)

where sub-matrix C _k1,k2 is an m × m covariance matrix between grids of the k ₁th and grids of the k ₂th level. When k ₁ = k ₂, the sub-matrix denotes the covariance of grids of a same level (diagonal element in (Eq. 8). K = 17 is the number of all levels.

The measurement operator H in (3) converts the model state variables to the observational variables. Its simplest form is a linear function to interpolate from the model grid to the observational stations, which is often used when the observed variable is the same as the analyzed state variable. In this study, the observed variable is SSTA and the analyzed state variable is STA at 17 levels including SSTA. For a full presentation, H = [H ₁ H ₂ …H _k H _k+1 … H _K ] is a L × M.matrix where H ₁ is a L × m matrix which interpolates model SSTA from model grids to SSTA on observation grids (L denoting the number of observation grids). H ₂, H ₃, …H _K are assumed, for simplicity, to be zero matrixes with the same dimension as matrix H ₁ Hence we can simplify terms in the Eq. 3 as follows:

$$ {\mathbf{P}}_e^f {\mathbf{H}}^T = \left[ {\begin{array}{*{20}c} {{\mathbf{C}}_{1,1} {\mathbf{H}}_1 } & {{\mathbf{C}}_{2,1} {\mathbf{H}}_1 } & \ldots & {{\mathbf{C}}_{K,1} {\mathbf{H}}_1 } \\ \end{array} } \right]^T $$

(9)

$$ {\mathbf{HP}}_e^f {\mathbf{H}}^T = \left[ {\begin{array}{*{20}c} {{\mathbf{H}}_1 {\mathbf{C}}_{1,1} } & {{\mathbf{H}}_1 {\mathbf{C}}_{1,2} } & \ldots & {{\mathbf{H}}_1 {\mathbf{C}}_{1,K} } \\ \end{array} } \right]{\mathbf{H}}^T = {\mathbf{H}}_1 {\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T $$

(10)

$$ ({\mathbf{d}} - {\mathbf{HX}}) = ({\mathbf{d}} - {\mathbf{H}}_1 {\mathbf{x}}_1 ) $$

(11)

Equation 3 could be rewritten as below:

$$ \left[ \begin{array}{c} {{\mathbf{x}}_1^a } \\ {{\mathbf{x}}_2^a } \\ \ldots \\ {\mathbf{x}}_k^a \\ \ldots \\ {{\mathbf{x}}_{16}^a } \\ {{\mathbf{x}}_{17}^a } \end{array} \right] = \left[ \begin{array}{c} {{\mathbf{x}}_1^f } \\ {{\mathbf{x}}_2^f } \\ \ldots \\ {\mathbf{x}}_k^f \\ \ldots \\ {{\mathbf{x}}_{16}^f } \\ {{\mathbf{x}}_{17}^f } \end{array} \right] + \left[ \begin{array}{c} {{\mathbf{C}}_{1,1} } \\ {{\mathbf{C}}_{2,1} } \\ \ldots \\ {\mathbf{C}}_{k,1} \\ \ldots \\ {{\mathbf{C}}_{16,1} } \\ {{\mathbf{C}}_{17,1} } \end{array} \right]{\mathbf{H}}_1 ^T ({\mathbf{H}}_1 {\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T + {\mathbf{R}}_e )^{ - 1} ({\mathbf{d}} - {\mathbf{H}}_1 {\mathbf{x}}_1^f ) $$

(12)

If C _1,1 is non-singular, the equation ( 12 ) could then be written:

$$ \left[ {\begin{array}{c} {{\mathbf{x}}_1^a } \\ {{\mathbf{x}}_2^a } \\ \ldots \\ {\mathbf{x}}_k^a \\ \ldots \\ {{\mathbf{x}}_{16}^a } \\ {{\mathbf{x}}_{17}^a } \\ \end{array} } \right] = \left[ {\begin{array}{c} {{\mathbf{x}}_1^f } \\ {{\mathbf{x}}_2^f } \\ \ldots \\ {\mathbf{x}}_k^f \\ \ldots \\ {{\mathbf{x}}_{16}^f } \\ {{\mathbf{x}}_{17}^f } \\ \end{array} } \right] + \left[ {\begin{array}{c} {\mathbf{I}} \\ {{\mathbf{B}}_{2,1} } \\ \ldots \\ {\mathbf{B}}_{k,1} \\ \ldots \\ {{\mathbf{B}}_{16,1} } \\ {{\mathbf{B}}_{17,1} } \\ \end{array} } \right]{\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T ({\mathbf{H}}_1 {\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T + {\mathbf{R}}_e )^{-1} ({\mathbf{d}}- {\mathbf{H}}_1 {\mathbf{x}}_1^f ) $$

(13)

Where I is an identity matrix,

$$ {\mathbf{B}}_{k,1} = {\mathbf{C}}_{k,1} {\mathbf{C}}_{1,1} ^{ - 1} $$

(14)

From (8), we have

$$ {\mathbf{x}}_1^a = {\mathbf{x}}_1^f + {\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T ({\mathbf{H}}_1 {\mathbf{C}}_{1,1} {\mathbf{H}}_1 ^T + {\mathbf{R}}_e )^{ - 1} ({\mathbf{d}} - {\mathbf{H}}_1 {\mathbf{x}}_1^f ) = {\mathbf{x}}_1^f + \Updelta {\mathbf{x}}_{\mathbf{1}} $$

(15)

$$ {\mathbf{x}}_k^a = {\mathbf{x}}_k^f + {\mathbf{B}}_{k,1} \Updelta {\mathbf{x}}_{\mathbf{1}} $$

(16)

Equation 16 shows that B _k,1 is a transfer matrix to project the corrections of SST (Δx ₁) to the subsurface of the kth level to correct its temperatures by B _k,1Δx ₁.

At each assimilation step, the transfer operator B _k,1 can be calculated by ensemble members with Eq. (8), thus it is a time-variant matrix. Theoretically, when the SST is perturbed with a random perturbation, B _k,1 could be obtained. However it was found that if we only perturbed SST, the C _k,1 would be very small due to a very small variation of the kth level temperature, leading to very small B _k,1. Thus in this study, the initial perturbations are exerted on not only SST but also the temperatures of other upper 4 model levels. From these initial perturbations, the ensemble is constructed. B _k,1 that is obtained by ensemble members of the prediction model, characterizes the physical relationship between SST and the subsurface temperature.

Usually the dimension of C _1,1 is far larger than the limited ensemble size, easily leading to a singular C _1,1. Theoretically it is still possible to compute the pseudo-inverse matrix of C _1,1. In practice, we directly use (8) to avoid computing the inverse matrix of C _1,1.