Evaluation of Ocean Data Assimilation in CAS-ESM-C: Constraining the SST Field

Climate change has been a hot topic in recent decades, both to the scientific community and to the public. Despite the considerable attention, it remains a debate as to the cause of the climate change in the 20th century: anthropogenic or internal variability of the climate system?

To solve this problem, it is not enough to use observational data alone, due to the short time duration of the instrumental record. Thus, coupled climate models (or earth system models, ESMs) are essential tools in the community to investigate related issues. To date, most scientists believe that anthropogenic greenhouse gases (GHGs) were very likely the main cause of the global warming since the second half of the 20th century. However, with respect to the so-called global warming hiatus at the beginning of the 21st century, the interdecadal Pacific Oscillation (IPO) or Pacific Decadal Oscillation (PDO) (Mantua et al., 1997; Mantua and Hare, 2002),which mostly reflect the internal variability of the air-sea coupled climate system on decadal or multi-decadal timescales, are regarded as the cause (e.g., Kosaka and Xie, 2013).

At regional scales, and taking the decadal variation of the East Asian summer monsoon (EASM) in the late 1970s as an example, the mechanism responsible for this transition is still disputable (Zhou et al., 2009). Two main viewpoints exist to explain the decadal variation of the EASM: the internal variability of the climate system (Yang and Lau, 2004; Li et al., 2010; Zhu et al., 2011; Qian and Zhou, 2014), or a response to external forcing (Menon et al., 2002; Kimoto, 2005; Jiang et al., 2013; Wang et al., 2013).

Although coupled models can respond to imposed external forcing (GHGs, aerosols, solar radiation, volcanos etc.), they cannot capture the internal variability of the climate system, such as the IPO. Thus, it is impossible for a CGCM to reproduce real climate fluctuations with actual timing without inputting observed information (Fujii et al., 2009). Previous studies show that the all-forcing runs of CMIP3 and CMIP5, which are the coupled historical simulations forced by all external forcing, cannot capture the decadal variation of the climate system, e.g., the EASM (Sun and Ding, 2008).

Owing to the fact that slow variations in a coupled atmosphere-ocean system tend to be controlled by the oceanic field, in this study, applying ocean data assimilation in a CGCM, we inputted the internal variability of climate system to provide a new scientific tool to the climate community, which can be used to investigate climate change and related issues. In this weakly coupled assimilation system, the air-sea coupling processes can be maintained and the actual climate internal fluctuations of the climate system (e.g., ENSO, PDO) can be inputted into the model. The main distinction of our work from previous studies on coupled data assimilation (e.g., Liu et al., 2013; Chang et al., 2013; Lu et al., 2015) is that our focus is mainly on whether ocean data assimilation in a coupled model can improve the climate simulation, e.g., the ENSO footprint and SST and precipitation relationship in the western North Pacific. In section 2, the model, data and assimilation scheme used are described. The climatology, ENSO global footprint and precipitation-SST relationship in the western North Pacific are then investigated in section 3, followed by a summary in section 4.

The model used in this study is a climate model version of the fully coupled ESM developed at the Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (CAS), named CAS-ESM-C (Sun et al., 2012). The atmospheric component (IAP AGCM4) in this coupled model is the latest version of the AGCM developed at the IAP. The horizontal resolution of IAP AGCM4 is 1.4^°× 1.4^°, with 26 vertical layers. The dynamical core of IAP AGCM4 is inherited from previous IAP models, while the physical processes are temporarily introduced from CAM, developed at NCAR. We are gradually substituting the physical packages in the model with our own schemes, such as the cloud-aerosol-radiation ensemble modeling system of (Liang and Zhang, 2013). In the model version used in this study, the Emanuel deep convection parameterization scheme was adopted (Emanuel, 1991). The performance of this model has been evaluated in simulating some basic features of the climate, e.g., the EASM (Dong et al., 2012; Zhang et al., 2013; Su et al., 2014; Dong and Xue, 2016).

The oceanic component (LICOM 1.0) is a global OGCM developed at the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG/IAP) (Liu et al., 2004). The land component (CLM3), sea ice component (CSIM5) and the coupler (coupler version 6) are all developed at NCAR (Dickinson et al., 2006; Sun et al., 2012; Dong et al., 2014). The performance of CAS-ESM-C in simulating ENSO has been systematically evaluated (Su et al., 2015). The model has also been used to investigate the decadal variation of the Aleutian low-Icelandic low relationship (Dong et al., 2014).

Three sets of observational gridded data [(1) to (3) below], one set of atmospheric reanalysis data (4) and one set of oceanic reanalysis data (5) were used in this study. The gridded datasets, regarded as observational data in the study, were produced from the objective analysis procedure and have their own bias compared with "real" observations.

(1) Monthly and daily mean GPCP data from 1979 to 2014, with a horizontal resolution of 2.5^°× 2.5^° (Huffman et al., 1997; Adler et al., 2003).

(2) Monthly mean SLP data provided by the Met Office Hadley Centre (HadSLP2), with a horizontal resolution of 5^°× 5^° (Allan and Ansell, 2006).

(3) The merged land air and SST anomaly analysis (er-ghcn-sst), which is based on data from the Global Historical Climatology Network of land temperatures and the ICOADS of SST data (Smith et al., 2008).

(4) The ERA-Interim dataset, which is an "interim" reanalysis to replace ERA-40 (Dee et al., 2011).

(5) SODA reanalysis data (Carton et al., 2000; Carton and Giese, 2008).

A set of quality controlled subsurface ocean temperature and salinity profiles and objective analyses (EN4) was also adopted to evaluate the performance of the assimilation system (Good et al., 2013).

Besides, the gridded daily OISST data (Reynolds et al., 2007) provided by NOAA, which spans from 1981 to 2014, were used to assimilate into the ocean model. The OISST dataset is an analysis constructed by combining observations from different platforms (e.g., satellites, ships, buoys) on a regular global grid. A spatially complete SST map is produced by interpolating to fill in gaps.

Figure 1.  Spatial distribution of SST (units: $^\circ$C) derived from the (a) observation, (b) CAS-ESM-C and (d) Assim, averaged from 1985 to 2014. (c) The SST differences between the base model and observation. (e) The SST differences between Assim and the observation. The observational data, assimilated into the model, are derived from AVHRR.

An ensemble optimal interpolation (EnOI) scheme was adopted to assimilate the OISST data into the ocean component of the coupled model (Evensen, 2003; Oke et al., 2008). The same assimilation scheme has been adopted in the establishment of a global ocean data assimilation system (Lu et al., 2014) and an ocean data assimilation system in the Indian Ocean and western Pacific Ocean (Yan et al., 2015). The description of the method provided here, in the following paragraph, reflects that of (Lu et al., 2014) and (Yan et al., 2015).

The following equations are solved to compute the analysis, based on the background field: \begin{eqnarray} {\omega}_{a}&=&{\omega}_{b}+{K}({\omega}_{o}-{H}{\omega}_{b}) ,(1)\\ {K}&=&({\rho}_{o}{P}){H}[{H}({\rho}_{o}{P}){H}^{T}+{R}]^{-1} ,(2)\\ {\omega}&=&[h_0,T,S,u,v]^{T} ,(3)\\ {P}&=&\alpha{A}'{A}'^{T}/(n-1) , (4)\end{eqnarray} where the subscripts a, b and o denote the analysis, background and observation, respectively, and the superscript T denotes the transpose of matrices. ω represents the state vector, including sea surface height (h₀), temperature (T), salinity (S), and zonal and meridional current (u and v). R denotes the observation error covariance matrix. H is the observation operator that interpolates from the model space to the observation space. ρ denotes a localizing correlation function used to remove the effects of sampling error due to the ensemble size being smaller than the dimension of the model space (Gaspari and Cohn, 1999). The circle between ρ and P denotes a Schur product. P represents the background error covariance matrix, which can be calculated from a matrix A', consisting of n ensemble members taken from the model state anomalies. As described in (Evensen, 2003), the stationary ensemble of model states adopted in EnOI is usually sampled during a long-term control run to estimate the structure of the background error covariance. Thus, before our assimilation experiment, a 200-yr control integration was conducted to derive the stationary ensemble of model states. The ensemble size n was 108 in this study. Furthermore, as in (Lu et al., 2014) and (Yan et al., 2015), the ensembles used in the assimilation were dependent on different months, in order to adequately describe the distinct characteristics of the oceanic current in different months. In Equation (4), α is a scalar parameter, which can alter the relative importance between the background error covariance matrix and the observation error covariance matrix.

In our assimilation experiment, OISST data were assimilated into the ocean model globally once every 7 days during the coupled control integration using the EnOI scheme described above. The time period was from 1981 to 2014. To assess the results of our assimilation experiment (hereafter Assim), the observational long-term mean SST field is shown in Fig. 1a, corresponding to the ω _o in Eq. (1). The background, corresponding to ω _b in Eq. (1), is shown in Fig. 1d. Model cold bias is located mainly in the tropics and midlatitudes, while warm bias is located mainly in the SH (Fig. 1e). The analysis field is shown in Fig.1b, corresponding to ω _a in Eq. (1). It is clear that, after the assimilation, much of the deviation in the SST field from the observation has been reduced, compared with that before the assimilation (Figs. 1b and c). The time series of the area-averaged forecast and analysis field in the NH, SH and tropical eastern Pacific (5^°S-5^°N, 150^°E-90^°W), compared with the observation, are further shown in Fig. 2. The forecast field can be pulled closer to the observation after the assimilation step, both on the global and regional scale.

During the assimilation step, although only SST observations are assimilated into the model, the subsurface sea surface height (h₀), temperature (T), salinity (S), and zonal (u) and meridional (v) current in each grid point are accordingly adjusted based on the background error covariance matrix [P in Eq. (2)]. The background error covariance s are based on the statistics of 108 static samples derived from a control integration of the coupled model. Figure 3 shows the climatological vertical profile of the global ocean temperature and salinity difference between the simulation (analysis and forecast) and the observation. For the vertical temperature profile, the model deviation has been reduced after the assimilation procedure, especially the upper layer (above 400 m) (Fig. 3a). And for the vertical salinity profile, the model deviation has also been reduced, in the layers around 50 m and 200 m (Fig. 3b). However, the model deviation has been increased (approximately 0.04 psu) after the assimilation procedure in the surface layer. Note that we only assimilated SST in this study. The degraded signal of the surface salinity field may be associated with the assimilation scheme in that the relationship between temperature and salinity is simply derived from static samples of the coupled model control run. In the future, we intend to assimilate more ocean fields (e.g., temperature and salinity profiles) in the coupled model to investigate whether this increased deviation in the surface can be further reduced. In the following section, we report the results of an evaluation of the assimilation system, including the climatology, ENSO global footprint, and precipitation-SST relationship in the western North Pacific.

Figure 2.  Area-averaged SST (units: $^\circ$C) in each assimilation step: (a) NH; (b) SH; (c) tropical eastern Pacific (5$^\circ$S-5$^\circ$N, 150$^\circ$E-90$^\circ$W). Red line denotes the observation; green line denotes the forecast; blue line denotes the analysis. The observed SSTs are derived from AVHRR data.

Figure 3.  Climatological vertical profile of the (a) global ocean temperature and (b) salinity difference between the simulation and observation, averaged from 1985 to 2014. Red (blue) line denotes the analysis (forecast) field. The observation is derived from a set of quality controlled subsurface ocean temperature and salinity profiles and objective analyses (EN4) (Good et al., 2013).

Figure 4.  The long-term averaged (1985-2014) subsurface temperature (units: $^\circ$C) distribution in the tropical Pacific (5$^\circ$S-5$^\circ$N, 120$^\circ$-280$^\circ$E) derived from (a) SODA, (b) CAS-ESM-C and (d) Assim, averaged from 1985 to 2014. (c) The differences between the base model and SODA. (e) The differences between Assim and SODA. The black line in (a, b, d) denotes the 20$^\circ$C isothermal line. The black boxes in (c, e) denote the bias in Assim is reduced compared to CAS-ESM-C.

Figure 4 shows the distribution of climatological tropical subsurface temperature in the tropical Pacific (5^°S-5^°N, 120^°-280^°E), from which the simulated structure of the thermocline in the tropical Pacific can be clearly seen. The results show that in the Assim experiment the evident cold bias in the base model, especially in the upper 100 m, can be largely reduced (black box in Figs. 4c and e). The same conclusion can also be drawn from the vertical temperature profile in Fig. 3. Meanwhile, note that in the Assim experiment the cold bias in the ocean model below 100 m still exists, and the simulated thermocline is still shallower than observed. In the future we intend to assimilate more fields (e.g., temperature and salinity profiles) besides SST in the model, to see whether this cold bias can be further reduced.

Figure 5.  Climatological spatial distribution of boreal summer (JJA) mean precipitation (units: mm d$^-1$) during 1985-2014 in the (a) observation and (b) assimilation, and (c) their differences. The observed precipitation is derived from monthly GPCP data. The pattern correlation coefficients (PCC) showed in (b) denotes the correlation coefficients between assimilation and observation.

Figure 6.  Climatological spatial pattern of (a, c) boreal summer (JJA) and (b, d) winter (DJF) mean SLP (units: hPa) in the (a, b) observation and (c, d) simulation, averaged from 1985 to 2014. The differences between Assim and the observation are shown in (e, f) for (e) boreal summer (JJA) and (f) boreal winter (DJF). The observed SLP is derived from monthly HadSLP2 data. The PCCs shown in (c, d) denote the correlation between the assimilation and observation.

Figure 5 shows the climatological spatial pattern of the June-August (JJA) mean precipitation in both the observation and simulation. The climatology is calculated from 1985 to 2014. The model can reasonably reproduce the JJA mean precipitation on the global scale, with a pattern correlation coefficient (PCC) of 0.72. The main rainfall center over the Asian-Australian monsoon region and the tropical eastern Pacific in the observation are reflects in the results of the Assim simulation. A wet bias is located mainly in the low latitudes of the South Pacific Ocean, which may lead to the double ITCZ problem——a common deviation that exists in state-of-the-art coupled models (Li and Xie, 2014; Zhang et al., 2015). Dry biases in the tropics and midlatitudes are still evident, especially in the northern midlatitudes (Fig. 5c).

The atmospheric general circulation is an important metric to evaluate a climate model. Figure 6 shows the simulated SLP in boreal summer and winter, and compares them with the observation. The Assim simulation can reasonably capture the main circulation systems, e.g., the North Pacific/Atlantic subtropical high in boreal summer, the continental low (high) pressure in Eurasia in boreal summer (winter), and SH subtropical high pressure (the Siberian high, Aleutian low, Icelandic low) in boreal summer (winter). The PCCs between the Assim simulation and the observation in boreal summer and winter are all approximately 0.95. The simulated climatological SLP is lower than observed in most areas in both boreal summer and winter (Figs. 6e and f). In the following subsection, we evaluate the performance of our system in simulating the variability of the climate system.

ENSO is the most important interannual variability of the climate system, which can leave "its footprint" around the globe. During the different phases of ENSO (developing, mature and decaying stages), it can influence the global climate very differently. To reasonably reproduce the ENSO global footprint is essential for climate models. In this subsection, we report the regression of surface air temperature (SAT), precipitation and the geopotential height in 850-hPa (Z850) fields on the observed Niño3.4 (5^°S-5^°N, 170^°-120^°W) index to investigate the ENSO global footprint in both the observation and simulation. Due to the fact that observed SST with actual timing is assimilated into our system, the observed Niño3.4 index can thus be used for regression-based diagnosis in both the observation and simulation.

Figure 7 shows the ENSO-related global SAT in El Niño developing autumn, mature winter, decaying spring and decaying summer. In the observation (Figs. 7a-d), a positive temperature anomaly in the eastern equatorial Pacific is associated with the evolution of El Niño, firstly developing before the mature phase in boreal winter and then decaying gradually. Positive temperature anomalies are also located in the North Indian Ocean, which may result from the teleconnection between the Pacific and Indian oceans through anomalous Walker circulation. A negative temperature anomaly exists mainly in the tropical western Pacific and Northwest/Southwest Pacific and North Atlantic in the midlatitudes. With respect to the continent, a west-east (northwest-southeast) dipole pattern is located in the Eurasia (North America) in El Niño developing autumn. In El Niño mature winter, a positive anomaly covers the most western part and eastern part of Eurasia and the northern part of North America. During the El Niño decaying phase (spring and summer), the eastern part of Eurasia and the North Pacific are mainly covered by a negative temperature anomaly, while positive anomalies are located in the northern part of North America.

The simulation reproduces the association between temperature anomalies and ENSO evolution, especially the positive anomaly in the tropical central and eastern Pacific and the negative anomaly in the western Pacific and North/South Pacific in the midlatitudes (Figs. 7e-h). Note that the decaying of El Niño is weaker in the simulation than in the observation; that is, the positive anomaly can exist until the following summer, with little weakening of the magnitude. This may be associated with deficiencies in the coupled model (Su et al., 2015).

Figure 7.  Regression of SAT (units: $^\circ$C) in El Niño (a, e) developing autumn, (b, f) mature winter, (c, g) decaying spring and (d, h) decaying summer, on the Niño3.4 index in the (a-d) observation and (e-h) simulation. The black box denotes the region in which the simulation is reasonable compared with the observation. The observational data are from the er-ghcn-sst dataset.

Figure 8.  As in Fig.7 but for precipitation (units: mm d$^-1$). The observations are monthly GPCP data.

For the continent, the model can capture the dipole temperature pattern in North America through the El Niño mature phase to the following spring (Figs. 7f-g). This may be a result of the fact that ENSO can influence North America through the Pacific-North America (PNA) teleconnection pattern, which can be reasonably reproduced in the model (figure not shown). Indeed, the simulation of large-scale circulation and teleconnection patterns is reasonable in most state-of-the-art climate models (e.g., Liang and Bradley, 2016). Besides, the model can reproduce the negative temperature anomaly in the eastern part of Eurasia and the North Pacific in El Niño decaying summer. The simulation of the temperature anomaly in Africa and South America is also quite reasonable.

Figure 8 shows the ENSO global footprint in precipitation. In the observation, the most evident influence is in the tropical region, including the positive response in the western Indian Ocean and central and eastern Pacific, and the negative response in the eastern Indian Ocean and western Pacific, throughout the evolution of El Niño (Figs. 8a-d). Note that the model can capture the main pattern of the precipitation response in the tropical area (Figs. 8e-h), except that the latitudinal width of the precipitation response is smaller than observed, which may be attributable to the deviation of the SST warming pattern from the observation during ENSO (Fig. 7). Besides, there are also notable biases in the model, e.g., the positive response in the central and eastern tropical Pacific extends much farther west than observed, and the precipitation response in ENSO decaying summer is stronger than that in the observation. Regionally, there are some exciting results insofar as the extratropical response is also reasonable, e.g., the model can reproduce the positive precipitation response in East Asia and the southern part of North America in ENSO mature winter and decaying spring. Thus, this assimilation system could be useful in seasonal and interannual climate prediction in some regions, e.g., East Asia.

Figure 9.  As in Fig.8 but for geopotential height at 850hPa (Z850) (units: 10 gpm). The observations are derived from the ERA-Interim dataset.

Figure 9 further shows the Z850 response to ENSO, which may reflect the anomalous general circulation associated with ENSO. In the observation (Figs. 9a-d), it can be seen that during ENSO evolution, the negative Z850 anomaly is mainly in the eastern tropical Pacific, and North/South Pacific in the midlatitudes, which is closely associated with the eastern Pacific surface warming and related atmospheric teleconnection pattern (e.g., the PNA/Pacific-South America teleconnection pattern). In the western Pacific and Maritime Continent region, there are positive Z850 anomalies. This west-high-east-low Z850 anomaly pattern is reminiscent of the negative phase of the Southern Oscillation associated with El Niño. The model can reasonably capture this anomalous Z850 pattern and its evolution with ENSO (Figs. 9e-h).

Figure 10.  (a) Observed and (b) simulated spatial pattern of the correlation coefficients between boreal summer (JJA) mean precipitation and SST anomalies. The observed precipitation is from GPCP, and the observed SST from AVHRR.

Air-sea coupling processes are important in the western North Pacific region and should be considered in numerical simulation (Wang et al., 2004, 2005; Zou and Zhou, 2013). However, AGCMs, which are driven by prescribed SST and sea ice, are usually adopted in short-term climate prediction (Bengtsson et al., 1993). Thus, the coupling between the atmosphere and ocean is neglected in such climate forecasts. One goal of our establishing of an assimilation system was to provide a tool, which considers air-sea coupling, to numerically predict the climate in the western North Pacific and East Asia area. Here, we briefly evaluate the relationship between summer precipitation and SST, which can reflect the air-sea coupled processes, especially in the western North Pacific region.

Figure 10 shows the correlation coefficients between boreal summer (JJA) precipitation anomalies and local SST anomalies. It is evident from the observation that anomalous summer precipitation in the central and eastern tropical Pacific is positively correlated with the local anomalous SST, indicating that it is the ocean surface temperature that drives the precipitation. Besides, positive correlation also exists in the subtropical Pacific in both the NH and SH. However, in the western North Pacific region, summer precipitation is negatively correlated with SST, indicating that the atmospheric feedback plays a major role in determining local SST. The assimilation system can reasonably reproduce this precipitation-SST relationship in the western North Pacific, confirming the results in (Wang et al., 2005) that air-sea coupling is rather important in this region.

Figure 11 further provides the lead-lag correlation between precipitation and SST in the western North Pacific region. Observation shows that the simultaneous correlation between precipitation and SST is nearly zero, while when SST leads (lags) precipitation by about 8-10 days, the correlation can reach 0.3 (-0.3), which reaches the 99% confidence level in terms of statistical significance. The correlation coefficient between the observation and simulation is 0.91. This relationship actually denotes the essential air-sea coupling in this region. Our assimilation system can reproduce this relationship, although the correlation coefficient is moderately larger than observed.

Figure 11.  Lead-lag correlations of 20-100-day band-pass filtered precipitation and SST in the western North Pacific region (5$^\circ$-35$^\circ$N, 110$^\circ$-155$^\circ$E). Red (blue) line denotes observation (simulation). The observed precipitation and SST are derived from daily GPCP data and AVHRR, respectively.

In this study, we developed a weaklycoupled data assimilation system, based on a fully coupled model (CAS-ESM-C), in which the observed SSTs are assimilated into the ocean model through an EnOI scheme. Several basic aspects derived from this system, including the climatology, ENSO global footprint and SST-rainfall relationship in the western North Pacific domain, have been assessed. Two notable advantages have been demonstrated: (1) The system can adequately reflect the atmospheric feedback as a fully coupled model; (2) The oceanic fields will not drift far away from the observation through ocean data assimilation. The conclusions can be summarized as follows:

(1) Assimilation module testing: The assimilation module can reduce the bias of both surface and subsurface temperature after assimilation. Compared to the observation, the spatial pattern and area-averaged time series of SST show much less bias in the analysis than the background. For the temperature profile, the model deviation has been reduced through the assimilation procedure, especially in the upper layer. However, degraded signals of surface salinity can still be found.

(2) The climatology: Our assimilation system can reasonably reproduce the JJA-mean precipitation on the global scale. The main rainfall centers in the observation were reproduced well by the Assim simulation experiment. However, the common double ITCZ problem still exists in this system. With respect to the atmospheric general circulation, the main circulation systems (e.g., the Siberian high, Aleutian low, and Icelandic low) are reasonably captured by our system.

(3) ENSO global footprint: The simulated SAT, precipitation and Z850 fields were regressed onto the observed Niño3.4 index to investigate the ENSO global footprint simulated by our assimilation system. The ENSO-related global SAT, Z850 and precipitation in El Niño developing autumn, mature winter, decaying spring and decaying summer were examined, and the main observed features were reasonably reproduced in the Assim simulation experiment.

(4) SST-rainfall relationship: In the western North Pacific region, air-sea coupling processes are crucial for precipitation simulation. The negative correlation between local summer SST and precipitation was reproduced in the Assim simulation experiment. Besides, the observed lead-lag correlations between intraseasonal precipitation and SST were also reflected well in the Assim simulation.