Importance of El Niño reproducibility for reconstructing historical CO₂ flux variations in the equatorial Pacific

Based on a set of climate simulations utilizing two kinds of Earth system models (ESMs) in which observed ocean hydrographic data are assimilated using exactly the same data assimilation procedure, we have clarified that the successful simulation of the observed air–sea CO₂ flux variations in the equatorial Pacific is tightly linked to the reproducibility of coupled physical air–sea processes. When an ESM with a weaker ENSO (El Niño–Southern Oscillations) amplitude than that of the observations was used for historical simulations with ocean data assimilation, the observed equatorial anticorrelated relationship between the sea surface temperature (SST) and the air–sea CO₂ flux on interannual to decadal timescales could not be represented. The simulated CO₂ flux anomalies were upward (downward) during El Niño (La Niña) periods in the equatorial Pacific. The reason for this was that the non-negligible correction term in the governing equation of ocean temperature, which was added via the ocean data assimilation procedure, caused an anomalous, spurious equatorial upwelling (downwelling) during El Niño (La Niña) periods, which brought more (less) subsurface layer water rich in dissolved inorganic carbon (DIC) to the surface layer. On the other hand, in the historical simulations where the observational data were assimilated into the other ESM with a more realistic ENSO representation, the correction term associated with the assimilation procedure remained small enough so as not to disturb an anomalous advection–diffusion balance for the equatorial ocean temperature. Consequently, spurious vertical transport of DIC and the resultant positively correlated SST and air–sea CO₂ flux variations did not occur. Thus, the reproducibility of the tropical air–sea CO₂ flux variability with data assimilation can be significantly attributed to the reproducibility of ENSO in an ESM. Our results suggest that, when using data assimilation to initialize ESMs for carbon cycle predictions, the reproducibility of the internal climate variations in the model itself is of great importance.

1 Introduction

Since the industrial revolution, vast quantities of greenhouse gases (e.g., CO₂) have been released into the atmosphere through human activities such as fossil fuel use and land use change. This increased atmospheric CO₂ concentration leads to global warming; however, both the oceanic and the terrestrial ecosystems absorb atmospheric CO₂ and are considered to work to relax the progress of the global warming (Sabine et al., 2004; Doney et al., 2009a, 2014; Le Quéré et al., 2009, 2010, 2016). Observation-based studies have reached the consensus that significant interannual variability in the air–sea CO₂ flux (hereafter CO2F) exists in some specific regions, such as the equatorial Pacific and the high latitudes of both hemispheres (e.g., Park et al., 2010; Valsala and Maksyutov, 2010; Landschützer et al., 2014; Rödenbeck et al., 2014), and the variation in CO2F associated with the El Niño–Southern Oscillation (ENSO) in the equatorial Pacific has been highlighted in many previous observation-based and simulation-based studies (Keeling and Revelle, 1985; Feely et al., 1997, 1999; Jones et al., 2001; Obata and Kitamura, 2003; McKinley et al., 2004; Patra et al., 2005). When an El Niño event occurs in the equatorial Pacific, the dissolved inorganic carbon (DIC) concentration in the surface layer decreases due to a reduction in the supply of cold, DIC-rich subsurface water to the surface layer compared with normal years, which stems from weaker equatorial upwelling associated with weaker trade winds (Le Borgne et al., 2002; Feely et al., 2004; Doney et al., 2009a, b). Correspondingly, the CO2F anomaly is downward during El Niño, and the inverse is seen during La Niña. Le Borgne et al. (2002) estimated that the upwelling of DIC-rich subsurface water accounts for up to 70 % of the CO2F variation in the equatorial Pacific, while the other 30 % is attributable to variation in the wind speed and biological processes. Accordingly, to estimate and predict variations in CO₂ uptake by the global ocean on timescales of several years, it would be informative to first consider the variations in the equatorial Pacific associated with ENSO.

The Paris Agreement is an agreement within the United Nations Framework Convention on Climate Change (UNFCCC, 2015) providing the framework of measures from 2021 to 2030 to act against climate change. The goal of the Paris Agreement is to restrict the rise in the global mean surface air temperature to well below 2 ^∘C relative to the preindustrial level. If greenhouse gas emissions continue to increase at their current rate, the Earth's surface will warm by 1.5 ^∘C within ∼ 20 years relative to the preindustrial state as reported in the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC, 2013). In this context, comprehensive understanding of the changes in the carbon cycle over previous years is essential for accurate predictions of the global carbon cycle, including natural variations, which will assist in evaluation of future CO₂ emission reductions (Kawamiya et al., 2020).

For future climate predictions, data assimilation procedures are incorporated into climate models in order to synchronize simulated climatic states in the model with observations, that is, the initialization of climate models. By incorporating data assimilation procedures into Earth system models (ESMs), it will be possible to reproduce and predict variations in biogeochemical properties (Brasseur et al., 2009; Tommasi et al., 2017a, b; Park et al., 2018). This includes an assessment of the predictability of CO2F on a decadal timescale for the global ocean (Li et al., 2016, 2019).

Focusing on CO2F fluctuations associated with ENSO in the equatorial Pacific, Dong et al. (2016) analyzed the results of the Earth system models (ESMs) that participated in the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012), which contributed to the AR5 of the Intergovernmental Panel on Climate Change (IPCC, 2013). They showed that only some ESMs could reproduce the observed anticorrelated relationship between SST and CO2F. This suggests that our understanding of ENSO and associated global carbon cycle variations are still insufficient. For reliable prediction of future CO₂ uptake on interannual to decadal timescales, it is necessary to understand coupled physical air–sea processes and the associated carbon cycle variations in the equatorial Pacific.

In this study, utilizing two kinds of ESMs in which observed ocean hydrographic data are assimilated, we attempted to identify the key processes to reproduce the observed historical air–sea CO₂ flux variations in the equatorial Pacific. The remainder of this paper is organized as follows: Sect. 2 provides a brief description of the models used in this study, the derived results are presented in Sect. 3, and a short discussion and summary are presented in Sect. 4.

2 Methods

2.1 Model description

In this study, we have conducted four experiments, NEW-assim, NEW, OLD-assim, and OLD. In NEW-assim and NEW, we used the MIROC-ES2L (Hajima et al., 2020), and in OLD-assim and OLD, we used the MIROC-ESM (Watanabe et al., 2011). The former is newly developed for CMIP Phase 6 (CMIP6; Eyring et al., 2016), whereas the latter is an official model of CMIP5. The physical core model of MIROC-ES2L is MIROC5.2, which is a minor update of MIROC5 (Watanabe et al., 2010; Tatebe et al., 2018), whereas the physical core model of MIROC-ESM is MIROC3m (K-1 model developers, 2004). The horizontal resolution of the atmospheric component of MIROC-ES2L (MIROC-ESM) has T42 spectral truncation (i.e., approximately 300 km) with 40 (80) vertical levels up to 3 hPa (0.003 hPa). The oceanic component of MIROC-ES2L has a horizontal tripolar coordinate system. In the spherical coordinate portion south of 63^∘ N, the longitudinal grid spacing is 1^∘, while the meridional grid spacing varies from approximately 0.5^∘ near the Equator to 1^∘ in midlatitude regions. There are 62 vertical levels in a hybrid σ–z coordinate system, the lowermost of which is located at a depth of 6300 m. The oceanic component of MIROC-ESM has a horizontal bipolar coordinate system: the longitudinal grid spacing of the oceanic component is approximately 1.4^∘, while the latitudinal grid intervals vary gradually from 0.5^∘ at the Equator to 1.7^∘ near both poles. There are 44 vertical levels in a hybrid σ–z coordinate system, the lowermost of which is located at a depth of 5300 m. The resolutions in MIROC-ES2L are higher than in MIROC-ESM. In particular, 31 (21) of the 62 (44) vertical layers in MIROC-ES2L (MIROC-ESM) are within the upper 500 m of depth. The increased number of vertical layers in MIROC-ES2L has been adopted in order to better represent the equatorial thermocline.

In NEW-assim and OLD-assim, we used the ESMs that incorporated the same simple scheme for ocean data assimilation, which comprised an incremental analysis update (IAU; Bloom et al., 1996; Huang et al., 2002). This technique is relatively simple compared with more elaborate techniques such as the ensemble Kalman filter and four-dimensional variational method, but it is widely used for decadal climate predictions (e.g., Mochizuki et al., 2010; Tatebe et al., 2012). A benefit of an IAU is its relatively low computational cost, which enables decadal- to centennial-scale integration and a variety of parameter sensitivity experiments. In an IAU, during the analysis interval from t=0 to t=τ, the governing equation including a correction term for temperature and salinity (X) that is written as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=\mathrm{adv}.+\mathrm{diff}.+F+{\displaystyle \frac{\mathit{\alpha }}{\mathit{\tau }}}\mathrm{\Delta }{X}^{\mathrm{a}},\end{array}$$

where adv. is the advection term; diff. is the diffusion term; F is the surface flux term; and the final term on the right-hand side is the correction term, where α is a constant, and ΔX^a is the analysis increment. We employed the values of τ=1 d and α=0.025, and the IAU was applied at depths between the sea surface and 3000 m (Tatebe et al., 2012). The analysis increment is calculated from ΔX^a=X^a(0) −X(0), where X^a(0) is the analysis and X(0) is the model first guess at t=0; this term is kept unchanged during the analysis interval from t=0 to t=τ. For X^a(0), we used observed anomalies with respect to the observed monthly mean climatology during the 1961–2000 period. For X(0), simulated anomalies in NEW-assim (OLD-assim) with respect to monthly mean climatology in NEW (OLD) were used. Such a scheme, often called “anomaly assimilation” or “anomaly initialization”, is also used in many previous studies (e.g., Smith et al., 2007; Keenlyside et al., 2008; Pohlmann et al., 2009; Li et al., 2016, 2019; Sospedra-Alfonso and Boer, 2020). The monthly objective analysis of ocean temperature and salinity (Ishii and Kimoto, 2009) is assimilated into the model as X^a, with linear interpolation to daily data. Because observed DIC concentrations are sparse in space and time, only ocean hydrographic data are used for data assimilation in the present study. Moreover, any atmospheric observations or reanalyses are not applied.

Both NEW and OLD are the exactly same as the historical simulations designated by the CMIP6 and CMIP5 protocols, respectively, with three ensemble members for each that are bifurcated from arbitrary years of the corresponding preindustrial control simulations. The ocean data assimilation experiments, NEW-assim and OLD-assim, are bifurcated from NEW and OLD at the year 1946, respectively, and they are integrated up to the year 2005. Note that the data assimilation experiments are driven with the same external forcings as in the historical simulations. In the later sections, the model results for 1961–2005 are analyzed.

2.2 Estimating pCO₂ change at the sea surface

CO2F depends on the difference in the CO₂ partial pressure between the sea and the air, i.e.,

$$\begin{array}{}\text{(2)}& \mathrm{CO}\mathrm{2}\mathrm{F}=K(p{\mathrm{CO}}_{\mathrm{2}}-p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{air}})(\mathrm{1}-\mathit{\gamma }),\end{array}$$

where pCO₂ (pCO${}_{\mathrm{2}}^{\mathrm{air}})$ is the CO₂ partial pressure in the sea (air); γ is the fraction of sea ice; and K=kα is the CO₂ gas transfer coefficient, where k represents the CO₂ gas transfer velocity (Wanninkhof, 1992, 2014), and α represents the solubility of CO₂ in seawater (Weiss, 1974). The CO₂ gas transfer velocity k is a function of wind speed and the Schmidt number (Wanninkhof, 1992). This study investigated the reproducibility of the anticorrelated relationship between CO2F and SST; therefore, the direction of the flux is important. As K does not affect the direction and the flux variation due to ENSO has larger amplitude in terms of pCO₂ than pCO${}_{\mathrm{2}}^{\mathrm{air}}$ (Dong et al., 2017), the direction of the flux is governed by the variation in pCO₂. Consequently, we evaluated the pCO₂ change at the sea surface in the equatorial Pacific.

Seawater pCO₂ values depend on temperature (T), salinity (S), the DIC concentration, and the total alkalinity (Alk); therefore, the change in pCO₂ can be expanded as follows:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}=\mathrm{C}\left(T\right)+\mathrm{C}\left(S\right)+\mathrm{C}\left(\mathrm{DIC}\right)+\mathrm{C}\left(\mathrm{Alk}\right)+\mathrm{Res}.,\end{array}$$

where C(X) = (∂pCO${}_{\mathrm{2}}/\partial X$)ΔX (X=T, S, DIC, Alk) is the pCO₂ change due to the change in X (X=T, S, DIC, Alk), and Res., which includes second-order terms (Takahashi et al., 1993), was estimated so that the left-hand and right-hand sides in Eq. (3) are equal in this study. In Sect. 3, we evaluate the CO2F and pCO₂ variations in the equatorial Pacific in NEW-assim, NEW, OLD-assim, and OLD, and we calculate each term in Eq. (3) for each experiment.

2.3 Observation and reanalysis dataset

To assess CO2F, ocean temperature, and wind speed of the model output, we used observation or reanalysis datasets. We used SOM-FFN as the CO2F dataset (Landschützer et al., 2016, 2017, 2018). It is an estimate based on the ocean surface CO₂ observation data collection, SOCATv3 (Bakker et al., 2016), and provides monthly data since 1982. It shows significant interannual variation in CO2F in some specific regions, such as the equatorial Pacific and the high latitudes of both hemispheres (Fig. S1). In Sect. 3, we focus on the CO2F in the Niño3 region (5^∘ S–5^∘ N, 150–90^∘ W) which shows notable variation in CO2F in the equatorial Pacific. This region is also the region of maximum variability for SST (Gill, 1980). The observational COBE-SST2 was used as the SST dataset (Ishii et al., 2005; Hirahara et al., 2014). The JRA-55 reanalysis (Kobayashi et al., 2015) was used as the wind speed dataset.

Figure 1 (a, c) Maps showing the correlation coefficient between monthly CO2F anomalies derived from SOM-FFN and that of (a) NEW-assim and (c) OLD-assim. The analysis period is from 1982 to 2005. The solid boxes show the Niño3 region (5^∘ S–5^∘ N, 90^∘–150^∘ W) and the dotted boxes show the Niño4 region (5^∘ S–5^∘ N, 160^∘ E–150^∘ W). (b, d) Time series of the detrended NINO3-SST (blue line) and NINO3-CO2F (red line, positive upward) anomalies simulated with (b) NEW-assim and (d) OLD-assim. Values plotted are the 1-year running mean, and the shading in panels (b) and (d) shows the ensemble spread (1σ). R denotes the correlation coefficients between the detrended ensemble mean NINO3-SST and NINO3-CO2F anomalies, with a 1-year running mean filter applied.

3 Results

3.1 CO₂ flux and pCO₂ anomalies in the Niño3 region

Horizontal maps of the correlation coefficients between simulated and observed CO2F values are shown in Fig. 1. The model output data were the ensemble mean and were linearly interpolated into the SOM-FFN grid. Note that the data were not detrended, and a 1-year running mean filter is applied to the monthly COF2 anomalies in the 1982–2005 period before calculating the correlation coefficients in accordance with the period for which the SOM-FFN dataset is available. CO2F in NEW-assim shows a positive correlation with SOM-FFN in the equatorial Pacific region (Fig. 1a) where significant interannual variations in CO2F are found (Fig. S1). On the other hand, CO2F in OLD-assim (Fig. 1c) is negatively correlated in the equatorial Pacific. The time series in the Niño3 region of both the 1-year running mean SST (hereafter, NINO3-SST) and CO2F (hereafter, NINO3-CO2F) anomalies simulated with NEW-assim (OLD-assim) are shown in Fig. 1b (Fig. 1d). Here, the data were detrended and monthly anomalies were calculated with respect to the 1971–2000 monthly mean climatology. The correlation coefficients between NINO3-SST and NINO3-CO2F anomalies in NEW-assim, OLD-assim, and the observations are −0.50, 0.44, and −0.75, respectively (Table 1). The results in NEW-assim are consistent with the observations, whereas those in OLD-assim are not. The correlation coefficients between the NINO3-SST and NINO3-CO2F anomalies in NEW and OLD are −0.85 and −0.67, respectively (Table 1 and Fig. S2). Note that OLD could capture the observed anticorrelated relationship between the NINO3-SST and NINO3-CO2F anomalies, but OLD-assim could not reproduce this relationship.

Table 1 Correlation coefficients between the detrended 1-year running mean NINO3-SST and NINO3-CO2F anomalies in NEW-assim, NEW, OLD-assim, OLD, and the observations. The correlations coefficients in NEW-assim, NEW, OLD-assim, and OLD are for the period from 1961 to 2005 (Figs. 1 and S2), and the correlations coefficient in the observations is for the period from 1982 to 2005.

As the vertical direction of CO2F is determined mainly by pCO₂ at the sea surface (see Eq. 2), we further estimated each term in Eq. (3) for each model output (Fig. 2). ∂pCO${}_{\mathrm{2}}/\partial X$ in the C(X) term in Eq. (3) (X=T, S, DIC, or Alk) was estimated based on the climatological annual mean T, S, DIC, and Alk at the sea surface within the Niño3 region in each experiment. ΔX in C(X) (ΔpCO₂ on the left-hand side of Eq. 3) is the variation in X (pCO₂) associated with ENSO and was calculated by averaging the monthly mean X (pCO₂) anomalies regressed on the NINO3-SST anomalies over the entire Niño3 region. Note that the NINO3-SST anomalies are standardized by the standard deviation. In the following, we describe the anomalies during El Niño periods, whereas the opposite applies during La Niña periods. In NEW-assim, NEW, and OLD, pCO₂ decreases because the effect of the decrease in pCO₂ with decreasing DIC concentrations is larger than that of the increase in pCO₂ with warming (Fig. 2). In OLD-assim, however, the effect of the increase in pCO₂ with warming is larger than that of OLD, and the decrease in pCO₂ with decreasing DIC concentrations is smaller than that of OLD, resulting in an increase in pCO₂. As noted in Sect. 1, previous studies (Le Borgne et al., 2002; Feely et al., 2004; Doney et al., 2009a, b) have shown that the variability in the upwelling during ENSO events dominates the equatorial Pacific CO2F variations through its regulation of DIC. In the following, we discuss the temperature and vertical velocity changes associated with ENSO along the Equator.

Figure 2 Changes in pCO₂ regressed onto the standardized NINO3-SST anomalies (ΔpCO₂) (dots) and its decomposition with changes in X (X=T, S, DIC, Alk), C(X), and Res. (Eq. 3) evaluated in NEW-assim, NEW, OLD-assim, and OLD. See the text for details regarding the calculation of ΔpCO₂, C(X), and Res.

3.2 DIC and vertical velocity changes

A cross section of the monthly ocean temperature anomalies regressed onto the standardized monthly mean NINO3-SST anomalies along the equatorial Pacific is presented in Figs. 3 and S3 in addition to the climatological annual mean depths of the 18, 20, and 22 ^∘C isotherms. Here, monthly temperature anomalies were calculated with respect to the 1971–2000 monthly mean climatology. The observational temperature anomalies and the climatological isotherms are derived from the monthly objective analysis of ocean temperature (Ishii and Kimoto, 2009). Amplitudes of the positive (negative) equatorial temperature anomalies in the upper (lower) layer of the eastern (western) equatorial Pacific in NEW are larger than in OLD and are closer to the observations. The intensity of ENSO, defined as the standard deviation of detrended 1-year running mean NINO3-SST anomalies from 1961 to 2005, is shown in Table 2. The intensity of ENSO in NEW is estimated to be 1.17 ^∘C (Table 2), which is a bit stronger than the observed value (0.80 ^∘C). On the other hand, the intensity of ENSO in OLD is 0.43 ^∘C, which is about half as large as the observed value. In addition, the climatological mean thermocline in NEW is tighter than in OLD and is closer to the observations. The improvement in ENSO reproducibility in NEW is mainly attributed to two updates in the model configuration. The first is the implementation of an updated plume model for cumulus convection with multiple cloud types where the lateral entrainment rate varies vertically depending on the surrounding environment (Chikira and Sugiyama, 2010). The state-dependent lateral entrainment affects the strength of the convectively induced coupled air–sea processes in the eastern tropical Pacific and, thus, the ENSO amplitude in the model. More details are described in Watanabe et al. (2010). The second is the reduction of numerical diffusion due to the introduction of a highly accurate tracer advection scheme in the ocean and an increase in the vertical resolution (Prather, 1986). The equatorial thermocline in the climatic-mean state of the tropical Pacific is more diffuse in OLD than in the observations, which is partly due to numerical diffusion, especially in vertical advection (Tatebe and Hasumi, 2010), and this model bias is much alleviated in NEW. Correspondingly, the so-called “thermocline mode” (e.g., Imada and Kimoto, 2006) becomes more effective and the ENSO amplitude becomes larger in NEW. As the ENSO amplitude in NEW is larger than in OLD, the variation in the equatorial trade winds, which causes anomalous equatorial vertical velocity, is also larger in NEW.

Table 2 The intensity and period of ENSO in NEW, OLD, and the observations calculated from the 1-year running mean NINO3-SST anomalies for the period from 1961 to 2005.

To assess the variations in zonal wind associated with ENSO, we estimated the 10 m zonal wind anomalies over the NINO4 region (5^∘ S–5^∘ N, 160^∘ E–150^∘ W; the dotted boxes in Fig. 1a) which are regressed onto the NINO3-SST anomalies (Table 3). The Niño4 region is the region of maximum variability for the equatorial trade winds (Fig. S4). Hereafter, the abovementioned regression coefficient is referred to as wind feedback (Guilyardi et al., 2006). The positive value of the wind feedback in NEW (0.92 m s⁻¹ K⁻¹) indicates westerly wind anomalies during El Niño, and this is consistent with the observational dataset, i.e., 1.02 m s⁻¹ K⁻¹. The wind feedback in OLD (0.46 m s⁻¹ K⁻¹) is about half of that in NEW and the observations, respectively.

Table 3 The wind feedback and the vertical velocity feedback in NEW-assim, NEW, OLD, and OLD-assim. The wind feedback is computed as the monthly 10 m zonal wind anomalies in the Niño4 region regressed onto the monthly NINO3-SST anomalies, and the vertical velocity feedback is the monthly vertical velocity anomalies at the depth of the 20 ^∘C isotherm in the Niño3 region regressed onto the monthly NINO3-SST anomalies. The wind feedback is also evaluated from the observation dataset.

NA: not available.

Figure 3 Anomalies of equatorial ocean temperature regressed onto the standardized NINO3-SST anomalies for NEW (a), OLD (b), and the observations (c). The contour interval is 0.1 ^∘C. The thick solid lines indicate the climatological-mean of the 18, 20, and 22 ^∘C isotherms.

Cross sections of the monthly upward water velocity and DIC concentration anomalies along the Equator regressed onto the standardized NINO3-SST anomalies in NEW (OLD) are shown in Fig. 4a and c (Fig. 4b and d), respectively. By reproducing wind feedback that is consistent with the observations, the westerly wind anomalies during El Niño periods in NEW (Fig. S4c) are comparable to that of the JRA-55 reanalysis (Fig. S4i), leading to a weakening of the upward vertical velocity of approximately $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ (Fig. 4a). This weakening of the upward vertical velocity causes a decrease in the surface DIC in the eastern equatorial Pacific during El Niño periods (Fig. 4c). In OLD, the smaller wind feedback and the associated smaller westerly wind anomalies than in the JRA-55 reanalysis (Fig. S4g) lead to a weakening of the upward vertical velocity of just 10⁻⁶ m s⁻¹ in the equatorial Pacific (Fig. 4b). Although the ENSO signal in OLD is weaker than in the observations, due to a decrease in the upward vertical velocity compared with normal years, the surface DIC concentration decreases during El Niño periods (Fig. 4d). This is consistent with Dong et al. (2016) and shows that OLD is able to qualitatively reproduce the negative correlation between the SST and DIC concentration anomalies in the eastern equatorial Pacific (Fig. S2b).

Figure 4 Anomalies of the equatorial vertical velocity (a, b) and DIC (c, d) regressed onto the standardized NINO3-SST anomalies for NEW (a, c) and OLD (b, d). Contour intervals are 0.5 × 10⁻⁶ m s⁻¹ in panels (a) and (b) and 2 μmol L⁻¹ in panels (c) and (d). The thick solid lines indicate the climatological mean of the 18, 20, and 22 ^∘C isotherms.

Figure 5 Equatorial temperature analysis increments (a, b), vertical velocity anomalies (c, d), and DIC anomalies (e, f) regressed onto the standardized NINO3-SST anomalies for NEW-assim (a, c, e) and OLD-assim (b, d, f). The contour intervals are 0.02 × 10^−6 ∘C s⁻¹ in panels (a) and (b), 0.5 × 10⁻⁶ m s⁻¹ in panels (c) and (d), and 2 µmol L⁻¹ in panels (e) and (f). The thick solid lines indicate the climatological mean of the 18, 20, and 22 ^∘C isotherms.

Next, we examined the correction term in temperature due to the data assimilation, i.e., the temperature analysis increment; the final term on the right-hand side of Eq. (1); and the variations in vertical velocity and DIC concentration. Anomalies of the monthly mean temperature analysis increments, the vertical velocity, and the DIC concentration along the Equator regressed onto the standardized NINO3-SST anomalies are shown in Fig. 5. The maximum absolute value of the equatorial temperature analysis increment in NEW-assim is found at depths of 10–40 m in the eastern equatorial Pacific, which is shallower than the depth of the thermocline (Fig. 5a). In NEW-assim, the wind feedback is 0.92 m s⁻¹ K⁻¹ (Table 3), which is of the same magnitude as that in NEW (0.92 m s⁻¹ K⁻¹), and the surface wind anomalies still show a similar pattern to that in NEW (Fig. S4a–d). The westerly wind anomalies in NEW-assim lead to a weakening of the upward vertical velocity along the Equator during El Niño periods (Fig. 5c). To assess the variation in the equatorial vertical velocity associated with ENSO, we estimated the anomalies of the vertical velocity at the depth of the 20 ^∘C isotherm (the depth of the thermocline) in the Niño3 region which are then regressed onto the NINO3-SST anomalies. Hereafter, the regression coefficient is referred to as the vertical velocity feedback. The vertical velocity feedback in NEW-assim is estimated to be $-\mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹ K⁻¹, which is not significantly different from that in NEW (−3.9 × 10⁻⁷ m s⁻¹ K⁻¹) (Table 3). The negative value of the vertical velocity feedback in NEW-assim indicates a weakening of the upward vertical velocity at the depth of the thermocline during El Niño periods in the eastern equatorial Pacific (Fig. 5c). The weakening of the upward vertical velocity causes a decrease in the supply of DIC-rich subsurface water to the surface layer, leading to a reduction in the surface DIC concentration (Fig. 5e). In OLD, the temperature variations associated with ENSO at the depth of the thermocline in the eastern equatorial Pacific are smaller than observed (see Fig. 3b and c), so that the correction term forces a 0.16 × 10^−6 ∘C s⁻¹ increase in the equatorial water temperature during El Niño periods in order to realize the observed temperature variations (Figs. 5b, S3b). The wind feedback in OLD-assim is 0.48 m s⁻¹ K⁻¹ (Table 3), which is the same as in OLD, and the map of the wind speed anomalies shows a similar pattern to that of OLD (Fig. S4e–h); however, the warming due to the data assimilation procedure during El Niño periods reduces the density, leading to low-pressure anomalies. This results in anomalous cyclonic circulation and convergence and, thus, an enhancement of the upward vertical velocity at the depth of the thermocline (Fig. 5d). The vertical velocity feedback in OLD-assim is 4.1 × 10⁻⁷ m s⁻¹ K⁻¹, which has an opposite sign to that in OLD, −4.9 × 10⁻⁷ m s⁻¹ K⁻¹ (Table 3). A positive value of the vertical velocity feedback indicates the enhancement of the upward vertical velocity at the depth of the thermocline during El Niño periods, which is inconsistent with the observations. This spurious enhancement of the upward vertical velocity during El Niño periods causes an increase in the surface DIC concentration (Fig. 5f), leading to a positive correlation between the SST and CO2F (Fig. 1d), contrary to observations. We have to note here that the vertical velocity distribution (Fig. 5c) is even different from NEW in NEW-assim due to the temperature analysis increment (Fig. 4a). As already discussed, the intensity of ENSO in NEW is slightly stronger than observed (Table 2). In addition, the period of ENSO, which is defined as the peak of the power spectrum of the 1-year running mean NINO3-SST, is 5.0 years in NEW, which is longer than the 3.5 years in the observations (see Table 2). Because the ENSO characteristics in NEW are not perfectly consistent with observations, the model nature, namely the responses of the vertical velocity and DIC concentration in ENSO, are still distorted by the temperature analysis increment, even in NEW-assim. This indicates that further model improvements are needed.

4 Discussion and summary

In the present study, comparing the results of two ESMs in which observed ocean hydrographic data are assimilated, we have clarified that the representation of the processes in the equatorial climate system is important to reproduce the observed anticorrelated relationship between the SST and CO2F in the equatorial Pacific. When the ocean temperature and salinity observations were assimilated into an ESM with weaker ENSO amplitude than the observations, the correction term in the governing equation of the ocean temperature, which was introduced in the data assimilation procedure, caused spurious upwelling (downwelling) anomalies along the Equator during El Niño (La Niña) periods, leading to an increased (decreased) supply of DIC-rich subsurface water to the surface layer. Due to the resultant increase (decrease) in the surface DIC concentration, the upward (downward) CO2F anomalies during El Niño (La Niña) periods were induced, which was inconsistent with observation. When the ocean temperature and salinity observations were assimilated into the other ESM with a rather realistic ENSO representation, the anticorrelated relationship between the SST and CO2F was reproduced.

Focusing on the CO2F fluctuations associated with ENSO in the equatorial Pacific, Dong et al. (2016) analyzed the results of the CMIP5 ESMs. They showed that only a portion of CMIP5 ESMs (including MIROC-ESM) could reproduce the observed anticorrelated relationship between the SST and CO2F. Bellenger et al. (2014) evaluated the reproducibility of ENSO in the CMIP5 models. They reported that most CMIP5 climate models and ESMs underestimate the amplitude of the wind stress feedback by 20 %–50 % and that only 20 % of CMIP5 models have a relative error within 25 % of the observed value. There are many ESMs where the ENSO characteristics and/or the SST–CO2F relationships are inconsistent with observations. Causes of this discrepancy should be addresses in future studies using methods such as multi-model analysis; moreover, process-based uncertainty estimation will also be required in initialized climate and carbon predictions as well as projections by ESMs.